# Clean up dataframe to keep only meaningful vars and drop NA, ready for LM
tidy_df <- function(src_df, ct, f, drop_islands=TRUE) {
clean_bei_df <- src_df %>%
select(CT_UID, all_of(all.vars(f))) %>%
drop_na() %>%
units::drop_units()
# Join geom from CT dataframe
df <- ct %>%
transmute(CT_UID = GeoUID) %>%
inner_join(clean_bei_df, by="CT_UID")
# Build neighborhood (with islands) and extract island IDs
ct_nb <- poly2nb(df)
if (drop_islands) {
islands <- lapply(ct_nb, min) %>% lapply(function(e) e == 0) %>% unlist %>% which()
# Cleanup subset DF to drop islands and NA
df <- df %>%
filter(!row_number() %in% islands) %>%
mutate(ct_no = row_number()) # Add CT number, matching row number
# Recompute neighborhood
ct_nb <- poly2nb(df)
# Clean up final dataframe (no island, no NA's, pure dataframe)
clean_df <- df %>%
as.data.frame() %>%
select(CT_UID, ct_no, all_of(all.vars(f)))
}
nbw = nb2listw(ct_nb, zero.policy = !drop_islands)
nbmatx = nb2mat(ct_nb, style = "B", zero.policy = !drop_islands)
list(df=clean_df, nbmatx=nbmatx, nbw=nbw, island_dropped=drop_islands)
}
This document builds on BEI_equity_paper.Rmd doc, which started as an analysis for the equity paper and moved to an exploration of the best modeling framework (simple OLS, spatial auto-regressive linear regression, spatial random effect model and INLA) to assess the relation between interventions and equity metrics. The conclusions of this exploration led to the following decisions:
spaMM package to model association; INLA could be the focus of another methodological paper which would compare both approachesThis document contains the analyses to be included in the BEI and equity paper.
We focus on obj #1: are urban interventions tend to be located in low SES neighborhoods. Here, the urban interventions considered are bike lanes and canopy/tree coverage and low SES ~ high Pampalon deprivation index. In a second step, we will look at the variations of UI and SES and their association.
Data extraction and pre-analyses for the paper looking at BEI and equity. BEI comprise bike lanes and canopy changes while equity is measured through Pampalon deprivation index. (Partially adapted from original work on BEI bike lanes – see bike_lane_stats.R and ReadMe.md.)
Paper is available here.
General processing steps:
In a second phase, these BEI changes will be linked to Pampalon index for 2016 (and 2011 ?)
UPDATE 2021-12-02 Following discussion with @Yan, add normalized bike line changes:
UPDATE 2021-12-08 Following discussion with @Yan
UPDATE 2022-01-14 Following discussion with @Ruben and @Yan
Year variable as continuous instead of category in LMEwSCOREMAT matches theoretical range of x-axis in graph (-.2 >> .2), might ponder the magnitude of the observed trendUPDATE 2022-02-10 Following gentrification meeting
wSCORESOC, which does not capture equity dimension (Meghan) [replaced by visible minority]UPDATE 2022-02-15
UPDATE 2022-06-28
UPDATE 2022-07-07
We use data categorized by Philippe Apparicio’s team who manually identified bike lanes for each census year since 1991. For this study, we limit ourselves to 2016, 2011 and 2006 census years.
On top of the original CT boundaries, three levels of buffer have been applied to the CT – 250m, 500m & 750m. Then the same series of processing steps (see above) have been applied to the buffers.
# Bike lanes, from Ph. Apparicio
reseau <- st_read(dsn="data/ReseauCyclableFinal.gdb", layer = "Reseau") # Already in NAD83 / MTM zone 8
## Reading layer `Reseau' from data source
## `/Users/benoit/WORKSPACE/gentrification_BEI_equity/data/ReseauCyclableFinal.gdb'
## using driver `OpenFileGDB'
## Simple feature collection with 82166 features and 72 fields
## Geometry type: GEOMETRY
## Dimension: XYZ, XYZM
## Bounding box: xmin: 266985.5 ymin: 5029251 xmax: 320986.1 ymax: 5062652
## z_range: zmin: 0 zmax: 43
## m_range: mmin: 0 mmax: 43
## Projected CRS: NAD83 / MTM zone 8
bike_lane <- reseau %>%
filter(An2016 == 1 | An2011 == 1) %>%
select(IdRte, ClsRte, Zone, starts_with("An"), starts_with("Typo_")) %>%
st_cast("MULTILINESTRING") # Get rid of a few MULTICURVE geometries
# CT boundaries for Montreal
CT16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='CT', geo_format = "sf") %>%
filter(Type == "CT") %>%
mutate(interact_aoi = CD_UID %in% c(2466, 2465, 2458)) %>% # Flag Montréal island, Laval and the South shore (Longueuil, St-Lambert, Brossard)
st_transform(st_crs(bike_lane))
## Reading geo data from local cache.
CT11 <- get_census(dataset='CA11', regions=list(CMA='24462'), level='CT', geo_format = "sf") %>%
filter(Type == "CT") %>%
st_transform(st_crs(bike_lane))
## Reading geo data from local cache.
compute_bikelane_by_area <- function(sf_areas, year_fld, typo_fld) {
# Compute length of bike lanes within each area.
# --
# Parameters:
# - sf_areas: sf class object defining the areas of interest, must have a GeoUID field
# - year_fld: field name specifying the year of interest, e.g. An2016
# - typo_fld: field name specifying the typology for the year of interest, e.g. Typo_2016
year_fld <- enquo(year_fld)
typo_fld <- enquo(typo_fld)
# Compute intersection of bike lanes with areas
bk <- bike_lane %>%
filter(!!year_fld == 1) %>%
st_intersection(sf_areas) %>%
mutate(bike_lane_length = st_length(.)) %>%
as.data.frame() %>%
group_by(GeoUID, !!typo_fld) %>%
summarise(bike_lane_length = sum(bike_lane_length)) %>%
ungroup() %>%
pivot_wider(names_from = !!typo_fld, names_prefix = "Bike_class", names_sort = TRUE,
values_from = bike_lane_length, values_fill = units::set_units(0, m)) %>%
mutate(Bike_lane_total = units::set_units(rowSums(select(., starts_with("Bike_class"))), m))
# Merge back into original sf_areas
bk <- sf_areas %>%
left_join(bk)
# Replace NA by 0, which occur in Bike_class length
bk[is.na(bk)] <- 0
return(bk)
}
compute_streetlength_by_area <- function(sf_areas) {
# Compute length of streets within each area.
# --
# Parameters:
# - sf_areas: sf class object defining the areas of interest, must have a GeoUID field
# Compute intersection of streets with areas
bk <- reseau %>%
st_cast("MULTILINESTRING") %>% # Get rid of a few MULTICURVE geometries
st_intersection(sf_areas) %>%
mutate(street_length = st_length(.)) %>%
as.data.frame() %>%
group_by(GeoUID) %>%
summarise(street_length = sum(street_length)) %>%
ungroup()
# Merge back into original sf_areas
bk <- sf_areas %>%
mutate(shape_area_km2 = units::set_units(st_area(.), 'km^2')) %>%
left_join(bk)
# Replace NA by 0
bk[is.na(bk)] <- 0
return(bk)
}
# Compute year 2016 and year 2011 bike lanes within 2016 CTs
# NB: contrary to the original work, we keep the same area of reference, i.e. 2016
bike_lane_by_CT16 <- compute_bikelane_by_area(CT16, An2016, Typo_2016)
bike_lane_by_CT11 <- compute_bikelane_by_area(CT16, An2011, Typo_2011)
bike_lane_by_CT06 <- compute_bikelane_by_area(CT16, An2006, Typo_2006)
# Compute total street length within CT/buffer
street_length_by_CT16 <- compute_streetlength_by_area(CT16)
# Reorganize data to have all data in one dataframe
bike_lane_changes <- CT16 %>%
left_join(select(as.data.frame(street_length_by_CT16), GeoUID, street_length, shape_area_km2), by="GeoUID") %>%
left_join(select(as.data.frame(bike_lane_by_CT16), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
left_join(select(as.data.frame(bike_lane_by_CT11), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2016ct", ".2011ct")) %>%
left_join(select(as.data.frame(bike_lane_by_CT06), GeoUID, starts_with("Bike_")), by="GeoUID")
# Compute ratio of bike lane vs street length
bike_lane_changes <- bike_lane_changes %>%
transmute(GeoUID = GeoUID,
interact_aoi = interact_aoi,
street_length = street_length,
street_length.hm = street_length / 100, # in hectometers, can be kept as float (not a dependent variable following a Poisson distribution)
Bike_lane_total.2016ct = Bike_lane_total.2016ct,
Bike_lane_total.2011ct = Bike_lane_total.2011ct,
Bike_lane_total.2006ct = Bike_lane_total,
Bike_lane_total.hm.2016ct = as.integer(round(Bike_lane_total.2016ct/100)),
Bike_lane_total.hm.2011ct = as.integer(round(Bike_lane_total.2011ct/100)),
Bike_lane_total.hm.2006ct = as.integer(round(Bike_lane_total/100)),
Bike_lane.by.street.2016ct = 100 * Bike_lane_total.2016ct / street_length, # In percents
Bike_lane.by.street.2011ct = 100 * Bike_lane_total.2011ct / street_length,
Bike_lane.by.street.2006ct = 100 * Bike_lane_total / street_length)
# Compute change between 2011 and 2016 (only for total bike lane length)
bike_lane_changes <- bike_lane_changes %>%
mutate(Bike_lane_diff.2011.2016ct = Bike_lane_total.2016ct - Bike_lane_total.2011ct)
# Normalize bike lane change by (i) street length and (ii) area
bike_lane_changes <- bike_lane_changes %>%
mutate(Bike_lane_diff.by.street.2011.2016ct = 100 * Bike_lane_diff.2011.2016ct / street_length)
# Save results
st_write(bike_lane_changes, dsn = "data/bike_length_changes.gpkg", delete_layer = TRUE)
## Deleting layer `bike_length_changes' using driver `GPKG'
## Writing layer `bike_length_changes' to data source
## `data/bike_length_changes.gpkg' using driver `GPKG'
## Writing 970 features with 15 fields and geometry type Multi Polygon.
# Clean up
rm(buf_CT16)
rm(bike_lane_by_CT11, bike_lane_by_CT16, bike_lane_by_CT06)
rm(street_length_by_CT16)
Check output for one specific dataset (Census tracts 2016, no buffer)
Bike lane length in 2016 within CTs, measured in meters
ggplot() +
geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_total.2016ct)), lwd=0) +
scale_fill_continuous(name = "Total length (m)")+
labs(title = "Length of bike lanes within 2016 CTs")
Absolute bike lane length change between 2011 and 2016, in meters
ggplot() +
geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.2011.2016ct)), lwd=0) +
scale_fill_gradient2(name = "Length (m)")+
labs(title = "Changes in bike lane between 2011 and 2016")
Relative bike lane length change between 2011 and 2016, normalized by street length within CT in 2016, expressed in %
\[Variation = \frac{Bike Lane_{2016}}{Street Length_{2016}} - \frac{Bike Lane_{2011}}{Street Length_{2016}}\]
ggplot() +
geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.by.street.2011.2016ct)), lwd=0) +
scale_fill_gradient2(name = "Variation (%)")+
labs(title = "Changes in bike lane between 2011 and 2016, normalized by street")
At the Census Tract level.
# CT level
ggplot(bike_lane_changes) +
geom_histogram(aes(as.numeric(Bike_lane_diff.2011.2016ct)), binwidth = 250) +
xlab("Difference of bike lane length between 2016 and 2011 | CT level")
summary(bike_lane_changes$Bike_lane_diff.2011.2016ct)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1377.7 0.0 0.0 262.1 189.7 14463.8
ggplot(bike_lane_changes) +
geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.2011.2016ct))) +
xlab("Difference of bike lane length between 2016 and 2011, normalized by street")
summary(bike_lane_changes$Bike_lane_diff.by.street.2011.2016ct)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -100.000 0.000 0.000 3.665 4.496 100.000 252
Canopy changes is based on data produced by CMM, using multispectral aerial imagery and lidar. In order to sync the observations with the census years, we focus on 2011 and 2019 with one extra observation point in 2015.
The processing steps are similar to the ones for the bike lanes:
# Codes du raster "espace vert"
# 0. No data (hors CMM)
# 1. NDVI < 0,3 et MNH < 3,0m = Minéral bas (route, stationnement, etc.)
# 2. NDVI < 0,3 et MNH ≥ 3,0m = Minéral haut (constructions)
# 3. NDVI ≥ 0,3 et MNH < 3,0m = Végétal bas (culture, gazon, etc.)
# 4. NDVI ≥ 0,3 et MNH ≥ 3,0m = Végétal haut (canopée)
# 5. Aquatique
# Load rasters into pg database for further processing
system('psql -d gentrif_bei -c "CREATE EXTENSION IF NOT EXISTS postgis"')
system('psql -d gentrif_bei -c "CREATE EXTENSION IF NOT EXISTS postgis_raster"')
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2019') IS NOT NULL;")) == 0) {
system("raster2pgsql -s 32188 -I -C -M data/canopy/2019/*.tif -F -t 1000x1000 canopee2019 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2019' already imported") }
## PG Raster 'canopee2019' already imported
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2017') IS NOT NULL;")) == 0) {
system("raster2pgsql -s 32188 -I -C -M data/canopy/2017/*.tif -F -t 1000x1000 canopee2017 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2017' already imported") }
## PG Raster 'canopee2017' already imported
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2015') IS NOT NULL;")) == 0) {
system("raster2pgsql -s 32188 -I -C -M data/canopy/2015/*.tif -F -t 1000x1000 canopee2015 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2015' already imported") }
## PG Raster 'canopee2015' already imported
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2011') IS NOT NULL;")) == 0) {
system("raster2pgsql -s 32188 -I -C -M data/canopy/2011/*.tif -F -t 1000x1000 canopee2011 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2011' already imported") }
## PG Raster 'canopee2011' already imported
# Resample to 10m as the original rasters have a 1m resolution, which is too high to allow for a swift processing
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2019_10m') IS NOT NULL;")) == 0) {
system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2019 mode=2\" -r mode -tr 10 10 data/canopy/canopee2019_10m.tif")
system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2019_10m.tif -F -t 100x100 canopee2019_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2019_10m' already imported") }
## PG Raster 'canopee2019_10m' already imported
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2017_10m') IS NOT NULL;")) == 0) {
system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2017 mode=2\" -r mode -tr 10 10 data/canopy/canopee2017_10m.tif")
system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2017_10m.tif -F -t 100x100 canopee2017_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2017_10m' already imported") }
## PG Raster 'canopee2017_10m' already imported
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2015_10m') IS NOT NULL;")) == 0) {
system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2015 mode=2\" -r mode -tr 10 10 data/canopy/canopee2015_10m.tif")
system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2015_10m.tif -F -t 100x100 canopee2015_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2015_10m' already imported") }
## PG Raster 'canopee2015_10m' already imported
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2011_10m') IS NOT NULL;")) == 0) {
system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2011 mode=2\" -r mode -tr 10 10 data/canopy/canopee2011_10m.tif")
system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2011_10m.tif -F -t 100x100 canopee2011_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2011_10m' already imported") }
## PG Raster 'canopee2011_10m' already imported
# Push CT16 to pg
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('ct16') IS NOT NULL;")) == 0) {
CT16 %>%
st_transform(crs = 32188) %>%
st_write(con_bei, "ct16",
layer_options = c("OVERWRITE=yes", "LAUNDER=true", "SPATIAL_INDEX=gist", "GEOMETRY_NAME=geom"))
system("psql -d gentrif_bei -c 'CREATE INDEX ON ct16 USING gist (geometry)'")
} else { message("PG Layer CT16 already imported") }
## PG Layer CT16 already imported
WITH cnt19 AS (
SELECT "GeoUID", "Population"
,(pvc).value, SUM((pvc).count) As total
FROM (SELECT "GeoUID", "Population"
,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
FROM canopee2019_10m
JOIN ct16 ON ST_Intersects(geometry, rast)
) As foo
GROUP BY "GeoUID", "Population", (pvc).value
),
canopee19 AS (
SELECT "GeoUID"
,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2019 -- area expressed in hectares
,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2019
,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2019
,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2019
FROM cnt19
WHERE value > 0 -- discard no data, including postgis raster no data
GROUP BY "GeoUID", "Population"
),
cnt17 AS (
SELECT "GeoUID", "Population"
,(pvc).value, SUM((pvc).count) As total
FROM (SELECT "GeoUID", "Population"
,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
FROM canopee2017_10m
JOIN ct16 ON ST_Intersects(geometry, rast)
) As foo
GROUP BY "GeoUID", "Population", (pvc).value
),
canopee17 AS (
SELECT "GeoUID"
,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2017
,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2017
,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2017
,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2017
FROM cnt17
WHERE value > 0 -- discard no data, including postgis raster no data
GROUP BY "GeoUID", "Population"
),
cnt15 AS (
SELECT "GeoUID", "Population"
,(pvc).value, SUM((pvc).count) As total
FROM (SELECT "GeoUID", "Population"
,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
FROM canopee2015_10m
JOIN ct16 ON ST_Intersects(geometry, rast)
) As foo
GROUP BY "GeoUID", "Population", (pvc).value
),
canopee15 AS (
SELECT "GeoUID"
,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2015
,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2015
,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2015
,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2015
FROM cnt15
WHERE value > 0 -- discard no data, including postgis raster no data
GROUP BY "GeoUID", "Population"
),
cnt11 AS (
SELECT "GeoUID", "Population"
,(pvc).value, SUM((pvc).count) As total
FROM (SELECT "GeoUID", "Population"
,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
FROM canopee2011_10m
JOIN ct16 ON ST_Intersects(geometry, rast)
) As foo
GROUP BY "GeoUID", "Population", (pvc).value
),
canopee11 AS (
SELECT "GeoUID"
,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2011
,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2011
,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2011
,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2011
FROM cnt11
WHERE value > 0 -- discard no data, including postgis raster no data
GROUP BY "GeoUID", "Population"
)
SELECT "GeoUID"
,round((st_area(geometry) / 10000)::numeric, 1) ct_area_h -- CT area in hectares
,COALESCE(area_esp_vert_2011, 0) area_esp_vert_2011
,coalesce(area_esp_vert_high_2011, 0) area_esp_vert_high_2011
,coalesce(pct_esp_vert_high_2011, 0) pct_esp_vert_high_2011
,coalesce(pct_esp_vert_2011, 0) pct_esp_vert_2011
,COALESCE(area_esp_vert_2015, 0) area_esp_vert_2015
,coalesce(area_esp_vert_high_2015, 0) area_esp_vert_high_2015
,coalesce(pct_esp_vert_high_2015, 0) pct_esp_vert_high_2015
,coalesce(pct_esp_vert_2015, 0) pct_esp_vert_2015
,COALESCE(area_esp_vert_2017, 0) area_esp_vert_2017
,coalesce(area_esp_vert_high_2017, 0) area_esp_vert_high_2017
,coalesce(pct_esp_vert_high_2017, 0) pct_esp_vert_high_2017
,coalesce(pct_esp_vert_2017, 0) pct_esp_vert_2017
,COALESCE(area_esp_vert_2019, 0) area_esp_vert_2019
,coalesce(area_esp_vert_high_2019, 0) area_esp_vert_high_2019
,coalesce(pct_esp_vert_high_2019, 0) pct_esp_vert_high_2019
,coalesce(pct_esp_vert_2019, 0) pct_esp_vert_2019
FROM ct16
FULL JOIN canopee19 USING ("GeoUID")
FULL JOIN canopee17 USING ("GeoUID")
FULL JOIN canopee15 USING ("GeoUID")
FULL JOIN canopee11 USING ("GeoUID");
Get it here
pampalon <- read.xlsx("data/Canada2016Pampalon/A-MSDIData_Can2016_eng/1. EquivalenceTableCanada2016_ENG.xlsx", sheet = 2) %>%
mutate(DA = as.character(DA)) %>%
select(DA, SCOREMAT, SCORESOC)
# 2016 DA boundaries for Montreal
DA16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='DA', geo_format = "sf") %>%
filter(Type == "DA") %>%
st_transform(st_crs(bike_lane))
## Reading geo data from local cache.
pampalon <- DA16 %>%
inner_join(pampalon, by = c("GeoUID" = "DA")) %>%
as.data.frame()
# Get Pampalon 2006
pampalon06 <- read.xlsx("data/Canada2006Pampalon/A-MSDIData_Can2006_eng/1. CorrespondenceTable_Can2006_eng.xlsx", sheet = 2) %>%
mutate(DA = as.character(DA)) %>%
select(DA, DAPOP2006, SCOREMAT, SCORESOC)
# Get LUT DA2006 <-> DA2011 from StatCan
lut_da.1 <- read.csv("data/2011_92-156_DA_AD_txt/2011_92-156_DA_AD.txt", colClasses = "character",
header = FALSE, col.names = c("DAUID2011.ADIDU2011", "DAUID2006.ADIDU2006", "DBUID2011", "DA_rel_flag")) %>%
select(!c(DBUID2011, DA_rel_flag)) %>%
unique()
# Link Pampalon 2011 to LUT and compute weighted mean of scores of Pampalon 2011
# NB: population numbers will diverge from reality when more than one DA is merged into one DA of next census
pampalon06.11 <- pampalon06 %>%
inner_join(lut_da.1, by = c("DA" = "DAUID2006.ADIDU2006")) %>%
group_by(DAUID2011.ADIDU2011) %>%
summarise(pop2006 = sum(DAPOP2006),
SCOREMAT.06 = weighted.mean(SCOREMAT, DAPOP2006, na.rm = TRUE),
SCORESOC.06 = weighted.mean(SCORESOC, DAPOP2006, na.rm = TRUE))
# Get Pampalon 2011
pampalon11 <- read.xlsx("data/Canada2011Pampalon/A-MSDIData_Can2011_eng/1. CorrespondenceTable_Can2011_eng.xlsx", sheet = 2) %>%
mutate(DA = as.character(DA)) %>%
select(DA, DAPOP2011, SCOREMAT, SCORESOC)
# Get LUT DA2011 <-> DA2016 from StatCan
lut_da <- read.csv("data/2016_92-156_DA_AD_csv/2016_92-156_DA_AD.csv", colClasses = "character") %>%
select(!c(DBUID2016.IDIDU2016, DA_rel_flag.AD_ind_rel)) %>%
unique()
# Link Pampalon 2011 to LUT, then to Pampalon 06 and finally compute weighted mean of scores of Pampalon 2011
pampalon11.16 <- pampalon11 %>%
inner_join(lut_da, by = c("DA" = "DAUID2011.ADIDU2011")) %>%
left_join(pampalon06.11, by =c("DA" = "DAUID2011.ADIDU2011")) %>%
group_by(DAUID2016.ADIDU2016) %>%
summarise(pop2011 = sum(DAPOP2011),
SCOREMAT = weighted.mean(SCOREMAT, DAPOP2011, na.rm = TRUE),
SCORESOC = weighted.mean(SCORESOC, DAPOP2011, na.rm = TRUE),
SCOREMAT.06 = weighted.mean(SCOREMAT.06, pop2006, na.rm = TRUE),
SCORESOC.06 = weighted.mean(SCORESOC.06, pop2006, na.rm = TRUE),
pop2006 = sum(pop2006))
# Then link Pampalon 2011 to 2016
pampalon <- pampalon %>%
left_join(pampalon11.16, by = c("GeoUID" = "DAUID2016.ADIDU2016"), suffix = c(".16", ".11"))
# Aggregate at the CT level
pampalon_CT <- pampalon %>%
group_by(CT_UID) %>%
summarise(wSCOREMAT.2016 = weighted.mean(SCOREMAT.16, Population, na.rm = TRUE),
wSCORESOC.2016 = weighted.mean(SCORESOC.16, Population, na.rm = TRUE),
wSCOREMAT.2011 = weighted.mean(SCOREMAT.11, pop2011, na.rm = TRUE),
wSCORESOC.2011 = weighted.mean(SCORESOC.11, pop2011, na.rm = TRUE),
wSCOREMAT.2006 = weighted.mean(SCOREMAT.06, pop2006, na.rm = TRUE),
wSCORESOC.2006 = weighted.mean(SCORESOC.06, pop2006, na.rm = TRUE))
# Clean up
rm(lut_da, lut_da.1, pampalon11.16, pampalon06.11, pampalon11, pampalon06)
# Display map
.pampalon_CT_geom <- CT16 %>%
left_join(pampalon_CT, by = c("GeoUID" = "CT_UID")) %>%
filter(interact_aoi)
.pampalon_data <- bi_class(.pampalon_CT_geom, x = wSCOREMAT.2016, y = wSCORESOC.2016, style = "quantile", dim = 3)
## Warning in classInt::classIntervals(bins_x, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes
## Warning in classInt::classIntervals(bins_y, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes
.map <- ggplot() +
geom_sf(data = .pampalon_data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
bi_scale_fill(pal = "DkBlue", dim = 3) +
labs(title = "Pampalon: material and social deprivation index") +
theme(panel.background = element_rect(fill = "white"),
#axis.ticks = element_blank(),
#axis.text = element_blank(),
panel.grid = element_line(color = "darkgray", size = 0.2))
.legend <- bi_legend(pal = "DkBlue",
dim = 3,
xlab = "Material ",
ylab = "Social ",
size = 8)
ggdraw() +
draw_plot(.map, 0, 0, 1, 1) +
draw_plot(.legend, 0.1, .7, 0.2, 0.2)
Using Ding metric computed on 5 year span.
# Load gentrified CTs, 5 year span (from repo gentrification_metrics)
ding <- list()
ding[["2016"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_16", quiet=TRUE) %>%
filter(cma_uid_16 == "24462") %>%
st_transform(st_crs(bike_lane))
ding[["2011"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_11", quiet=TRUE) %>%
filter(cma_uid_11 == "24462") %>%
st_transform(st_crs(bike_lane))
ding[["2006"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_06", quiet=TRUE) %>%
filter(cma_uid_06 == "24462") %>%
st_transform(st_crs(bike_lane))
.ding_map <- ding[["2016"]] %>%
left_join(select(as.data.frame(CT16), GeoUID, interact_aoi), by = c("ct_uid_16" = "GeoUID")) %>%
filter(interact_aoi)
ggplot(data = .ding_map) +
geom_sf(aes(fill = gentrified_2016_2011, colour=gentrifiable_2011)) +
scale_fill_manual(values = c("gray", "red", "darkgray"), name = "Gentrified in 2016") +
scale_colour_manual(values = c("darkgray", "darkred", "darkgray"), name = "Gentrifiable in 2011") +
labs(title = "Census tract gentrification status in 2016")
Introduced here as a proposition, nothing acted (2022-02-04)
# Visible Minority
# - v_CA16_3954: Total - Visible minority for the population in private households - 25% sample data (Total)
# - v_CA16_3957: Total visible minority population (Total)
# Low income (LIM-AT)
# - v_CA16_2540: Prevalence of low income based on the Low-income measure, after tax (LIM-AT) (%) (Total)
equity_ct16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='CT', geo_format = "sf",
vectors = c("v_CA16_3954", "v_CA16_3957", "v_CA16_2540")) %>%
filter(Type == "CT") %>%
transmute(CT_UID = GeoUID,
vis_minority_2016 = `v_CA16_3957: Total visible minority population` / `v_CA16_3954: Total - Visible minority for the population in private households - 25% sample data` * 100,
low_income_2016 = `v_CA16_2540: Prevalence of low income based on the Low-income measure, after tax (LIM-AT) (%)`)
## Reading vectors data from local cache.
## Reading geo data from local cache.
# Visible Minority
# - v_CA11N_457: CA 2011 NHS, Total population in private households by visible minority (Total)
# - v_CA11N_460: CA 2011 NHS, Total population in private households by visible minority, Total visible minority population (Total)
# Low income (LIM-AT)
# - v_CA11N_2606: CA 2011 NHS, Prevalence of low income in 2010 based on after-tax low-income measure % (Total)
equity_ct11 <- get_census(dataset='CA11', regions=list(CMA='24462'), level='CT', geo_format = "sf",
vectors = c("v_CA11N_457", "v_CA11N_460", "v_CA11N_2606")) %>%
filter(Type == "CT") %>%
transmute(CT_UID = GeoUID,
vis_minority_2011 = `v_CA11N_460: Total visible minority population` / `v_CA11N_457: Total population in private households by visible minority` * 100,
low_income_2011 = `v_CA11N_2606: Prevalence of low income in 2010 based on after-tax low-income measure %`)
## Reading vectors data from local cache.
## Reading geo data from local cache.
# Visible Minority
# - v_CA06_1302: Total population by visible minority groups
# - v_CA06_1303: Total population by visible minority groups, Total visible minority population
# Low income (LIM-AT)
# - v_TX2006_551: After-tax low income status of tax filers and dependents (census family low income measure, CFLIM-AT) for couple and lone parent families by family composition, 2006 | All family units | Persons in Low Income | % - Total
equity_ct06 <- get_census(dataset='CA06', regions=list(CMA='24462'), level='CT', geo_format = "sf",
vectors = c("v_CA06_1302", "v_CA06_1303", "v_TX2006_551")) %>%
filter(Type == "CT") %>%
transmute(CT_UID = GeoUID,
vis_minority_2006 = `v_CA06_1303: Total visible minority population` / `v_CA06_1302: Total population by visible minority groups - 20% sample data` * 100,
low_income_2006 = `v_TX2006_551: % - Total`)
## Reading vectors data from local cache.
## Reading geo data from local cache.
equity_ct <- st_join(equity_ct16, equity_ct11, left=TRUE, largest=TRUE, suffix=c("", "_2011")) %>% # join on largest overlap, to overcome mismatch in CT UID
st_join(equity_ct06, left=TRUE, largest=TRUE, suffix=c("", "_2006")) %>%
data.frame()
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# cleanup
rm(equity_ct11, equity_ct16, equity_ct06)
# Display map
.equity_CT_geom <- CT16 %>%
left_join(equity_ct, by = c("GeoUID" = "CT_UID")) %>%
filter(interact_aoi)
.equity_data <- bi_class(.equity_CT_geom, x = vis_minority_2016, y = low_income_2016, style = "quantile", dim = 3)
## Warning in classInt::classIntervals(bins_x, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes
## Warning in classInt::classIntervals(bins_y, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes
.map <- ggplot() +
geom_sf(data = .equity_data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
bi_scale_fill(pal = "Brown", dim = 3) +
labs(title = "Equity metrics: % of visible minority and % of low-income household") +
theme(panel.background = element_rect(fill = "white"),
#axis.ticks = element_blank(),
#axis.text = element_blank(),
panel.grid = element_line(color = "darkgray", size = 0.2))
.legend <- bi_legend(pal = "Brown",
dim = 3,
xlab = "Vis. Minority ",
ylab = "Low-Income ",
size = 8)
ggdraw() +
draw_plot(.map, 0, 0, 1, 1) +
draw_plot(.legend, 0.1, .7, 0.2, 0.2)
All variables + outcome linked at the CT level + quintiles of SES variables
.bike_lane_changes <- bike_lane_changes %>%
as.data.frame() %>%
select(GeoUID, starts_with("street_length"), ends_with("ct", ignore.case = FALSE)) %>%
select(GeoUID, starts_with("street_length"), starts_with("Bike_lane")) # Drop individual category lane length
bei_df <- CT16 %>%
as.data.frame() %>%
transmute(CT_UID = GeoUID,
CD_UID = CD_UID,
CSD_UID = CSD_UID,
interact_aoi = interact_aoi,
zone = case_when(CD_UID == "2466" ~ "Montreal",
CD_UID == "2458" ~ "Longueuil",
CD_UID == "2465" ~ "Laval",
TRUE ~ "Other"),
Population = Population) %>%
left_join(pampalon_CT, by="CT_UID") %>%
left_join(select(as.data.frame(ding$`2016`), ct_uid_16, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_16")) %>%
left_join(select(as.data.frame(ding$`2011`), ct_uid_11, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_11")) %>%
left_join(select(as.data.frame(ding$`2006`), ct_uid_06, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_06")) %>%
left_join(select(as.data.frame(equity_ct), !c("geometry", "CT_UID_2011", "CT_UID_2006")), by="CT_UID") %>%
left_join(.bike_lane_changes, by=c("CT_UID" = "GeoUID")) %>%
left_join(as.data.frame(rename_with(esp_vert_ct, ~ paste0(., "ct"))), by=c("CT_UID" = "GeoUIDct")) %>%
# Compute quintile of SES variables
mutate(wSCOREMAT.2006.Q = ntile(wSCOREMAT.2006, 5),
wSCOREMAT.2011.Q = ntile(wSCOREMAT.2011, 5),
wSCOREMAT.2016.Q = ntile(wSCOREMAT.2016, 5),
low_income_2006.Q = ntile(low_income_2006, 5),
low_income_2011.Q = ntile(low_income_2011, 5),
low_income_2016.Q = ntile(low_income_2016, 5),
vis_minority_2006.Q = ntile(vis_minority_2006, 5),
vis_minority_2011.Q = ntile(vis_minority_2011, 5),
vis_minority_2016.Q = ntile(vis_minority_2016, 5)) %>%
units::drop_units()
head(bei_df)
write.csv(bei_df, "data/_results/bei_equity.csv", na="", row.names = FALSE)
# keep only interact CT
# UPDATE 2022-07-07 ##########
# and recompute quintiles for the new AOI
bei_df_aoi <- filter(bei_df, interact_aoi) %>%
mutate(wSCOREMAT.2006.Q = ntile(wSCOREMAT.2006, 5),
wSCOREMAT.2011.Q = ntile(wSCOREMAT.2011, 5),
wSCOREMAT.2016.Q = ntile(wSCOREMAT.2016, 5),
low_income_2006.Q = ntile(low_income_2006, 5),
low_income_2011.Q = ntile(low_income_2011, 5),
low_income_2016.Q = ntile(low_income_2016, 5),
vis_minority_2006.Q = ntile(vis_minority_2006, 5),
vis_minority_2011.Q = ntile(vis_minority_2011, 5),
vis_minority_2016.Q = ntile(vis_minority_2016, 5))
Included variables:
CT_UID: 2016 Census Tract IDCD_UID: 2016 Census DivisionCSD_UID: 2016 Census Subdivisioninteract_aoi: Does CT belong to INTERACT study area?zone: Montreal / Laval / Longueuil, etc. / OtherPopulation: 2016 Population within CTct_area_m2ct: Area of CT, in square metersgentrified_2016_2011: Is the CT gentrified in 2016?gentrifiable_2011: Is the CT candidate to gentrification in 2011?gentrified_2011_2006: Is the CT gentrified in 2011gentrifiable_2006: Is the CT candidate to gentrification in 2006gentrified_2006_2001: Is the CT gentrified in 2006gentrifiable_2001: Is the CT candidate to gentrification in 2001wSCORESOC.2016: Social deprivation index in 2016 (population weighted)wSCOREMAT.2016: Material deprivation index in 2016 (population weighted)wSCOREMAT.2016.Q: Quintile of material deprivation index in 2016 (population weighted)wSCORESOC.2011: Social deprivation index in 2011 (population weighted)wSCOREMAT.2011: Material deprivation index in 2011 (population weighted)wSCOREMAT.2011.Q: Quintile of aterial deprivation index in 2011 (population weighted)wSCORESOC.2006: Social deprivation index in 2006 (population weighted)wSCOREMAT.2006: Material deprivation index in 2006 (population weighted)wSCOREMAT.2006.Q: Quintile of material deprivation index in 2006 (population weighted)vis_minority_2016: % of visible minority in CT 2016vis_minority_2016.Q: Quintile of % of visible minority in CT 2016low_income_2016: prevalence of low income in CT 2016low_income_2016.Q: Quintile of prevalence of low income in CT 2016vis_minority_2011: % of visible minority in CT 2011vis_minority_2011.Q: Quintile of % of visible minority in CT 2011low_income_2011: prevalence of low income in CT 2011low_income_2011.Q: Quintile of prevalence of low income in CT 2011vis_minority_2006: % of visible minority in CT 2006vis_minority_2006.Q: Quintile of % of visible minority in CT 2006low_income_2006: prevalence of low income in CT 2006low_income_2006.Q: Qintile of prevalence of low income in CT 2006Bike_lane_total.{2016|2011}ct: total length of bike lanes, in 2016 or 2011, within CTBike_lane_total.hm.{2016|2011}ct: total length of bike lanes, in hectometers (as integers), in 2016 or 2011, within CTBike_lane.by.street.{2016|2011}ct: % of bike lanes compared to streets, in 2016 or 2011, within CTBike_lane_diff.2011.2016ct: change in total length of bike lanes between 2011 and 2016, within CTBike_lane_diff.by.street.2011.2016ct: change in total length of bike lanes between 2011 and 2011, normalized by street length, within CTBike_lane_diff.by.area.2011.2016ct: change in total length of bike lanes between 2011 and 2011, normalized by area, within CTarea_esp_vert_{2011|2015|2019}ct: area of green space in 2011, 2015 or 2019 within CT (in hectares)area_esp_vert_high_{2011|2015|2019}ct: same as above, except for trees (high canopy)pct_esp_vert_{2011|2015|2019}ct: % of green space in 2011, 2015 or 2019 within CTpct_esp_vert_high_{2011|2015|2019}ct: same as above, except for trees (high canopy)pct_esp_vert_diff{2011|2015}.{2015|2019}ct: change in % of green space between 2011 and 2015, 2011 and 2019 as well as 2011 and 2019, within CTINTERACT study area ~ Montréal, Laval, Longueuil, Brossard, St-Lambert
.bei_df_long <- bei_df_aoi %>%
select(CT_UID, CD_UID, starts_with("wSCORE")) %>%
select(!ends_with('.Q')) %>%
pivot_longer(!c(CT_UID, CD_UID))
ggplot(.bei_df_long, aes(value)) +
geom_histogram() +
facet_wrap(~name) #, scales = "free")
## Warning: Removed 86 rows containing non-finite values (stat_bin).
.bei_df_long <- bei_df_aoi %>%
select(CT_UID, CD_UID, starts_with("vis_minority"), starts_with("low_income")) %>%
select(!ends_with('.Q')) %>%
pivot_longer(!c(CT_UID, CD_UID))
ggplot(.bei_df_long, aes(value)) +
geom_histogram() +
facet_wrap(~name) #, scales = "free")
## Warning: Removed 85 rows containing non-finite values (stat_bin).
.bei_df_long <- bei_df_aoi %>%
select(CT_UID, CD_UID, starts_with("gentrif")) %>%
select(!ends_with('.Q')) %>%
pivot_longer(!c(CT_UID, CD_UID))
ggplot(.bei_df_long, aes(value)) +
geom_bar() +
facet_wrap(~name) #, scales = "free", ncol = 3)
.bei_df_long <- bei_df_aoi %>%
filter(interact_aoi) %>%
select(CT_UID, CD_UID, matches("^Bike_lane_total.*ct$"), matches("^area_esp_vert.*ct$")) %>%
pivot_longer(!c(CT_UID, CD_UID))
ggplot(.bei_df_long, aes(value)) +
geom_histogram() +
facet_wrap(~name, scales = "free", ncol = 4)
.bei_df_long <- bei_df_aoi %>%
filter(interact_aoi) %>%
select(CT_UID, CD_UID, matches("^Bike_lane.by.*ct$"), matches("^pct_esp_vert.*ct$")) %>%
pivot_longer(!c(CT_UID, CD_UID))
ggplot(.bei_df_long, aes(value)) +
geom_histogram() +
facet_wrap(~name, scales = "free", ncol = 3)
Looking at objective #1: do urban interventions tend to be located in low SES neighborhoods?. We look at \[Urban Condition_{2011} = f(SES_{2011})\] as well as \[Urban Condition_{2011} = f(Gentrification_{2011 \to 2016})\]
Here \(UrbanCondition\) means the state of the urban environment features, such as length of bike lanes, greenness coverage, etc. at one specific moment. This needs to be distinguished from \(UrbanIntervention\), which accounts for the changes in the \(UrbanConditions\) between two years (see below).
We fit LMM models, with spatial random effect using spaMM package.
Bike lane ratio to streets at the Census tract level. We use an offset, as explained in Anderson’s presentation, to handle the denominator and model the Bike lane length in hectometers using a Poisson distribution.
Resorting to a hurdle model, to accommodate the zero-inflated distribution.
f <- Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2011ct != 0) ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.3714 0.36381 -6.518 7.111e-11
## wSCOREMAT.2011.Q -0.2359 0.09349 -2.523 1.164e-02
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1377511
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 1.489
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -354.34
ci.lmm.0 <- confint(res.lmm.zero, 'wSCOREMAT.2011.Q')
## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q
## -0.44869647 -0.04680899
res.lmm.0 <- update(res.lmm.zero, . ~ . -wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm.zero)
## chi2_LR df p_value
## p_v 6.028834 1 0.01407403
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
plot(simout, title = "DHARMa residual | Logit")
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2011ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.04308 0.09569 -21.351 0.00000
## wSCOREMAT.2011.Q -0.06073 0.02808 -2.163 0.03058
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1340074
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.3915
## # of obs: 469; # of groups: ct_no, 469
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1735.752
ci.lmm.cnt <- confint(res.lmm.cnt, 'wSCOREMAT.2011.Q')
## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q
## -0.116508140 -0.005344248
res.lmm.0 <- update(res.lmm.cnt, . ~ . -wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm.cnt)
## chi2_LR df p_value
## p_v 4.616437 1 0.03166695
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4188.184
Interpretation of the wSCOREMAT.2011.Q coefficient (significant in both sub-models): - we have \(e^{-0.236}\) = 0.79 (95%CI: 0.638 – 0.954), which means each increase of 1 quintile in the Material Score 2011 leads to a 21% decrease the chance of having a bicycle lane in the CT - we have \(e^{-0.061}\) = 0.941 (95%CI: 0.89 – 0.995), which means in CTs with bike lanes, each increase of 1 quintile in the Material Score 2011 leads to a 5.9% decrease the ratio of bike lanes to street in 2011.
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 at the Census tract level, using a Poisson distribution
f <- area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -0.59856 0.04009 -14.931 0
## wSCOREMAT.2011.Q -0.09327 0.01094 -8.527 0
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.138379
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.05515
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2410.162
ci.lmm <- confint(res.lmm, 'wSCOREMAT.2011.Q')
## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q
## -0.11472563 -0.07177056
# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 68.11096 1 1.110223e-16
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.093}\) = 0.911 (95%CI: 0.892 – 0.931), which means each increase of 1 quintile in the Material Score leads to a 8.9% decrease in the percentage of greenspace area within CT in 2011.
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 at the Census tract level, using a Poisson distribution
f <- area_esp_vert_high_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.2846 0.06041 -21.27 0.000e+00
## wSCOREMAT.2011.Q -0.1294 0.01655 -7.82 5.329e-15
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1376847
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1306
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2154.052
ci.lmm <- confint(res.lmm, 'wSCOREMAT.2011.Q')
## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q
## -0.1619265 -0.0968720
# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 58.1149 1 2.475797e-14
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{NA}\) = NA (95%CI: 0.851 – 0.908), which means each increase of 1 quintile in the Material Score leads to a 12.1% decrease in the percentage of tree canopy area within CT in 2011.
Bike lane ratio to streets (in %) at the Census tract level
f <- Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ vis_minority_2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2011ct != 0) ~ vis_minority_2011.Q + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.536 0.33404 -7.592 3.153e-14
## vis_minority_2011.Q -0.212 0.08987 -2.359 1.834e-02
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1377585
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 1.356
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -354.9056
ci.lmm.0 <- confint(res.lmm.zero, 'vis_minority_2011.Q')
## lower vis_minority_2011.Q upper vis_minority_2011.Q
## -0.40490772 -0.02914871
res.lmm.0 <- update(res.lmm.zero, . ~ . - vis_minority_2011.Q)
LRT(res.lmm.0, res.lmm.zero)
## chi2_LR df p_value
## p_v 5.144237 1 0.02332365
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
plot(simout, title = "DHARMa residual | Logit")
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2011ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.03046 0.09576 -21.204 0.00000
## vis_minority_2011.Q -0.06903 0.02867 -2.408 0.01605
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1297798
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.3876
## # of obs: 469; # of groups: ct_no, 469
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1730.437
ci.lmm.cnt <- confint(res.lmm.cnt, 'vis_minority_2011.Q')
## lower vis_minority_2011.Q upper vis_minority_2011.Q
## -0.12556369 -0.01248152
res.lmm.0 <- update(res.lmm.cnt, . ~ . - vis_minority_2011.Q)
LRT(res.lmm.0, res.lmm.cnt)
## chi2_LR df p_value
## p_v 5.704475 1 0.01692172
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4178.686
Interpretation of the vis_minority_2011.Q coefficient: - model zero-hurdle (non significant): we have \(e^{-0.212}\) = 0.809 (95%CI: 0.667 – 0.971), which means each increase of 1 quintile in the proportion of visible minority 2011 leads to a NA% decrease the chance of having a bicycle lane in the CT - model count: we have \(e^{-0.069}\) = 0.933 (95%CI: 0.882 – 0.988), which means in CTs with bike lanes, each increase of 1 quintile in the Material Score 2011 leads to a 6.7% decrease the ratio of bike lanes to street in 2011.
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -0.68172 0.04080 -16.711 0.00e+00
## vis_minority_2011.Q -0.06827 0.01151 -5.931 3.02e-09
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1385835
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.05926
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2421.774
ci.lmm <- confint(res.lmm, 'vis_minority_2011.Q')
## lower vis_minority_2011.Q upper vis_minority_2011.Q
## -0.09084917 -0.04560224
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.068}\) = 0.934 (95%CI: 0.913 – 0.955), which means each increase of 1 quintile in the proportion of visible minority leads to a 6.6% decrease in the percentage of greenspace area within CT in 2011.
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.3612 0.06080 -22.39 0.000e+00
## vis_minority_2011.Q -0.1094 0.01722 -6.35 2.152e-10
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380368
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1366
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2159.726
ci.lmm <- confint(res.lmm, 'vis_minority_2011.Q')
## lower vis_minority_2011.Q upper vis_minority_2011.Q
## -0.14321003 -0.07552651
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.109}\) = 0.896 (95%CI: 0.867 – 0.927), which means each increase of 1 quintile in the proportion of visible minority leads to a 10.4% decrease in the percentage of tree canopy area within CT in 2011.
Gentrified CT between 2011 and 2016
Bike lane ratio to streets (in %) at the Census tract level
f <- Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ gentrified_2016_2011 + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2011ct != 0) ~ gentrified_2016_2011 + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -3.12397 0.1952 -16.003 0.0000
## gentrified_2016_2011TRUE -0.06627 0.2473 -0.268 0.7887
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1378747
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 1.496
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -357.4465
ci.lmm.0 <- confint(res.lmm.zero, 'gentrified_2016_2011TRUE')
## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE
## -0.6062718 0.4522833
res.lmm.0 <- update(res.lmm.zero, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm.zero)
## chi2_LR df p_value
## p_v 0.06253228 1 0.8025374
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
plot(simout, title = "DHARMa residual | Logit")
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2011ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.26118 0.05223 -43.290 0.0000
## gentrified_2016_2011TRUE 0.09494 0.08002 1.186 0.2355
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1267498
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.3946
## # of obs: 469; # of groups: ct_no, 469
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1732.647
ci.lmm.cnt <- confint(res.lmm.cnt, 'gentrified_2016_2011TRUE')
## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE
## -0.06921123 0.25966646
res.lmm.0 <- update(res.lmm.cnt, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm.cnt)
## chi2_LR df p_value
## p_v 1.285916 1 0.2568019
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4188.186
Interpretation of the gentrified_2016_2011 coefficient: none significant - model zero-hurdle (non significant): we have \(e^{-0.066}\) = 0.936 (95%CI: 0.545 – 1.572) - model count: we have \(e^{0.095}\) = 1.1 (95%CI: 0.933 – 1.296)
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -0.8304 0.02578 -32.205 0.000e+00
## gentrified_2016_2011TRUE -0.1964 0.03343 -5.875 4.234e-09
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1385641
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.05997
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2421.862
ci.lmm <- confint(res.lmm, 'gentrified_2016_2011TRUE')
## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE
## -0.2620145 -0.1306189
# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 33.56416 1 6.895166e-09
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.196}\) = 0.822 (95%CI: 0.769 – 0.878), which means gentrified CTs exhibit a 17.8% decrease in the percentage of greenspace area within CT in 2011.
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.6372 0.03951 -41.442 0.000000
## gentrified_2016_2011TRUE -0.1505 0.04990 -3.015 0.002568
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380114
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1457
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2174.734
ci.lmm <- confint(res.lmm, 'gentrified_2016_2011TRUE')
## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE
## -0.24862544 -0.05240792
# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 9.025009 1 0.002663107
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.15}\) = 0.86 (95%CI: 0.78 – 0.949), which means gentrified CTs exhibit a 14% decrease in the percentage of tree canopy area within CT in 2011
UPDATE 2022-08-18: adding the same models, with discrete SES variables
f <- Bike_lane_total.hm.2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ factor(wSCOREMAT.2011.Q) + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2011ct != 0) ~ factor(wSCOREMAT.2011.Q) +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.6782 0.3458 -7.7452 9.548e-15
## factor(wSCOREMAT.2011.Q)2 -0.4105 0.3926 -1.0455 2.958e-01
## factor(wSCOREMAT.2011.Q)3 -0.1724 0.4091 -0.4215 6.734e-01
## factor(wSCOREMAT.2011.Q)4 -0.3350 0.4134 -0.8105 4.176e-01
## factor(wSCOREMAT.2011.Q)5 -1.1552 0.4092 -2.8232 4.754e-03
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1376613
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 1.49
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -352.3199
res.lmm.0 <- update(res.lmm.zero, . ~ . -factor(wSCOREMAT.2011.Q))
LRT(res.lmm.0, res.lmm.zero)
## chi2_LR df p_value
## p_v 10.06913 4 0.03927914
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
plot(simout, title = "DHARMa residual | Logit")
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2011ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.0448 0.08497 -24.064 0.000000
## factor(wSCOREMAT.2011.Q)2 -0.1713 0.10504 -1.631 0.102979
## factor(wSCOREMAT.2011.Q)3 -0.2997 0.10741 -2.790 0.005265
## factor(wSCOREMAT.2011.Q)4 -0.1538 0.11014 -1.397 0.162501
## factor(wSCOREMAT.2011.Q)5 -0.3097 0.12349 -2.508 0.012140
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1315181
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.3866
## # of obs: 469; # of groups: ct_no, 469
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1733.093
res.lmm.0 <- update(res.lmm.cnt, . ~ . -factor(wSCOREMAT.2011.Q))
LRT(res.lmm.0, res.lmm.cnt)
## chi2_LR df p_value
## p_v 9.933888 4 0.04155618
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4184.826
Transforming the wSCOREMAT.2011.Q coefficients (discrete):
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.069 |
| factor(wSCOREMAT.2011.Q)2 | 0.663 |
| factor(wSCOREMAT.2011.Q)3 | 0.842 |
| factor(wSCOREMAT.2011.Q)4 | 0.715 |
| factor(wSCOREMAT.2011.Q)5 | 0.315 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.129 |
| factor(wSCOREMAT.2011.Q)2 | 0.843 |
| factor(wSCOREMAT.2011.Q)3 | 0.741 |
| factor(wSCOREMAT.2011.Q)4 | 0.857 |
| factor(wSCOREMAT.2011.Q)5 | 0.734 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 at the Census tract level, using a Poisson distribution
f <- area_esp_vert_2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -0.71469 0.03647 -19.594 0.000e+00
## factor(wSCOREMAT.2011.Q)2 -0.08137 0.03997 -2.036 4.176e-02
## factor(wSCOREMAT.2011.Q)3 -0.12132 0.04238 -2.863 4.201e-03
## factor(wSCOREMAT.2011.Q)4 -0.22252 0.04335 -5.133 2.850e-07
## factor(wSCOREMAT.2011.Q)5 -0.40739 0.04796 -8.494 0.000e+00
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1383781
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.05402
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2406.448
# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - factor(wSCOREMAT.2011.Q))
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 75.5386 4 1.554312e-15
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Transforming the wSCOREMAT.2011.Q coefficients (discrete):
# Zero model
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.489 |
| factor(wSCOREMAT.2011.Q)2 | 0.922 |
| factor(wSCOREMAT.2011.Q)3 | 0.886 |
| factor(wSCOREMAT.2011.Q)4 | 0.801 |
| factor(wSCOREMAT.2011.Q)5 | 0.665 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 at the Census tract level, using a Poisson distribution
f <- area_esp_vert_high_2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2011ct ~ factor(wSCOREMAT.2011.Q) + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.4370 0.05487 -26.188 0.000e+00
## factor(wSCOREMAT.2011.Q)2 -0.1328 0.06019 -2.207 2.732e-02
## factor(wSCOREMAT.2011.Q)3 -0.1644 0.06378 -2.578 9.934e-03
## factor(wSCOREMAT.2011.Q)4 -0.3365 0.06576 -5.116 3.113e-07
## factor(wSCOREMAT.2011.Q)5 -0.5602 0.07293 -7.682 1.565e-14
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1377046
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1292
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2151.268
# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - factor(wSCOREMAT.2011.Q))
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 63.68247 4 4.874989e-13
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Transforming the wSCOREMAT.2011.Q coefficients (discrete):
# Zero model
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.238 |
| factor(wSCOREMAT.2011.Q)2 | 0.876 |
| factor(wSCOREMAT.2011.Q)3 | 0.848 |
| factor(wSCOREMAT.2011.Q)4 | 0.714 |
| factor(wSCOREMAT.2011.Q)5 | 0.571 |
Bike lane ratio to streets (in %) at the Census tract level
f <- Bike_lane_total.hm.2011ct ~ factor(vis_minority_2011.Q) + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ factor(vis_minority_2011.Q) + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2011ct != 0) ~ factor(vis_minority_2011.Q) +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.96762 0.3115 -9.526787 0.00000
## factor(vis_minority_2011.Q)2 0.24833 0.3889 0.638526 0.52313
## factor(vis_minority_2011.Q)3 -0.30809 0.3714 -0.829465 0.40684
## factor(vis_minority_2011.Q)4 0.00379 0.3925 0.009655 0.99230
## factor(vis_minority_2011.Q)5 -0.90644 0.3872 -2.340770 0.01924
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1376616
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 1.41
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -352.3024
res.lmm.0 <- update(res.lmm.zero, . ~ . - factor(vis_minority_2011.Q))
LRT(res.lmm.0, res.lmm.zero)
## chi2_LR df p_value
## p_v 10.35059 4 0.03491857
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
plot(simout, title = "DHARMa residual | Logit")
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2011ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ factor(vis_minority_2011.Q) + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2011ct ~ factor(vis_minority_2011.Q) + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.13756 0.08458 -25.274 0.00000
## factor(vis_minority_2011.Q)2 0.06116 0.10160 0.602 0.54720
## factor(vis_minority_2011.Q)3 -0.15828 0.10811 -1.464 0.14317
## factor(vis_minority_2011.Q)4 -0.19282 0.11075 -1.741 0.08168
## factor(vis_minority_2011.Q)5 -0.20731 0.12665 -1.637 0.10165
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1300438
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.3847
## # of obs: 469; # of groups: ct_no, 469
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1728.974
res.lmm.0 <- update(res.lmm.cnt, . ~ . - factor(vis_minority_2011.Q))
LRT(res.lmm.0, res.lmm.cnt)
## chi2_LR df p_value
## p_v 8.631541 4 0.07099879
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4176.552
Transforming the wSCOREMAT.2011.Q coefficients (discrete):
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.051 |
| factor(vis_minority_2011.Q)2 | 1.282 |
| factor(vis_minority_2011.Q)3 | 0.735 |
| factor(vis_minority_2011.Q)4 | 1.004 |
| factor(vis_minority_2011.Q)5 | 0.404 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.118 |
| factor(vis_minority_2011.Q)2 | 1.063 |
| factor(vis_minority_2011.Q)3 | 0.854 |
| factor(vis_minority_2011.Q)4 | 0.825 |
| factor(vis_minority_2011.Q)5 | 0.813 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2011ct ~ factor(vis_minority_2011.Q) + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2011ct ~ factor(vis_minority_2011.Q) + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2011ct ~ factor(vis_minority_2011.Q) + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -0.76472 0.03621 -21.118 0.000e+00
## factor(vis_minority_2011.Q)2 -0.05421 0.04138 -1.310 1.902e-01
## factor(vis_minority_2011.Q)3 -0.10559 0.04297 -2.457 1.401e-02
## factor(vis_minority_2011.Q)4 -0.14248 0.04466 -3.190 1.423e-03
## factor(vis_minority_2011.Q)5 -0.30573 0.04934 -6.196 5.793e-10
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1385852
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.05874
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2418.952
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Transforming the wSCOREMAT.2011.Q coefficients (discrete):
# Zero model
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.465 |
| factor(vis_minority_2011.Q)2 | 0.947 |
| factor(vis_minority_2011.Q)3 | 0.900 |
| factor(vis_minority_2011.Q)4 | 0.867 |
| factor(vis_minority_2011.Q)5 | 0.737 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2011ct ~ factor(vis_minority_2011.Q) + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2011ct ~ factor(vis_minority_2011.Q) + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2011ct ~ factor(vis_minority_2011.Q) + offset(log(ct_area_hct)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.48392 0.05407 -27.443 0.000e+00
## factor(vis_minority_2011.Q)2 -0.09687 0.06175 -1.569 1.167e-01
## factor(vis_minority_2011.Q)3 -0.17639 0.06414 -2.750 5.961e-03
## factor(vis_minority_2011.Q)4 -0.29448 0.06714 -4.386 1.154e-05
## factor(vis_minority_2011.Q)5 -0.45506 0.07409 -6.142 8.157e-10
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380636
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1366
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2159.165
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
Transforming the wSCOREMAT.2011.Q coefficients (discrete):
# Zero model
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.227 |
| factor(vis_minority_2011.Q)2 | 0.908 |
| factor(vis_minority_2011.Q)3 | 0.838 |
| factor(vis_minority_2011.Q)4 | 0.745 |
| factor(vis_minority_2011.Q)5 | 0.634 |
No new models here, as gentrified variable is already used as a pure categorical variable.
Looking at objective #1 | do urban interventions tend to be located in low SES neighborhoods?. We look at \[Urban Intervention_{2011 \to 2016} = f(SES_{2011})\] as well as \[Urban Intervention_{2011 \to 2016} = f(Gentrification_{2011 \to 2016})\]
Here \(Urban Intervention\) means the changes in the urban environment features, such as variation of bike lane length, greenness coverage, etc.
Apparently, change in rates are rather modeled as a pre-post model, as suggested in Jones’s answer. Hence, we model change in interventions between 2011 and 2016 as \[UC_{2016} = \beta_{0} + \beta_{1} * UC_{2011} + \beta_{2} * SES\]
Mostly redundant with multi-level approach below.
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -5.12247 0.8914 -5.7468 9.094e-09
## wSCOREMAT.2011.Q 0.05264 0.2153 0.2445 8.069e-01
## Bike_lane.by.street.2011ct 3.93126 2.5484 1.5426 1.229e-01
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct -0.39374 0.5959 -0.6607 5.088e-01
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1377292
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 5.737
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -149.3657
# NB: not converging | but not needed, use Wald p-value instead or LRT
# ci.lmm.0.1 <- confint(res.lmm.0, "wSCOREMAT.2011.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Bike_lane.by.street.2011ct")
# ci.lmm.0.3 <- confint(res.lmm.0, "wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct")
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
plot(simout)
}
# Plotting interaction || MEANINGLESS
# .t <- clean_bei$df
# .t["predict"] <- as.vector(predict(res.lmm.0))
#
# ggplot(.t, aes(y=predict, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
# geom_smooth(method = "lm", alpha=.1) +
# geom_point() +
# labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -2.6735086 0.085129 -31.40532
## wSCOREMAT.2011.Q -0.0084696 0.023916 -0.35414
## Bike_lane.by.street.2011ct 0.0537190 0.004433 12.11790
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct 0.0001078 0.001406 0.07667
## p-value
## (Intercept) 0.0000
## wSCOREMAT.2011.Q 0.7232
## Bike_lane.by.street.2011ct 0.0000
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct 0.9389
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1502892
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1125
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1884.428
# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "wSCOREMAT.2011.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Bike_lane.by.street.2011ct")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct")
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))
g_mat_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
geom_smooth(method = "lm", alpha=.1, linetype="dashed") +
labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_bk.hrdl
# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4079.588
Interpretation of the wSCOREMAT.2011.Q coefficient | Zero model (not significant): we have \(e^{0.053}\) = 1.054 Interpretation of the wSCOREMAT.2011.Q coefficient | Count model (not significant): we have \(e^{-0.008}\) = 0.992
Coefficient interpretation (marginal + interaction) as relative risk.
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ wSCOREMAT.2011.Q*pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.449918 0.0493086 -29.405 0.000e+00
## wSCOREMAT.2011.Q -0.061732 0.0149037 -4.142 3.442e-05
## pct_esp_vert_2011ct 0.015957 0.0008800 18.132 0.000e+00
## wSCOREMAT.2011.Q:pct_esp_vert_2011ct 0.001138 0.0003134 3.633 2.804e-04
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1364178
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004538
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2010.108
ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.2011.Q")
## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q
## -0.09102817 -0.03232795
ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_2011ct")
## lower pct_esp_vert_2011ct upper pct_esp_vert_2011ct
## 0.01421609 0.01772789
ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.2011.Q:pct_esp_vert_2011ct")
## lower wSCOREMAT.2011.Q:pct_esp_vert_2011ct
## 0.0005235603
## upper wSCOREMAT.2011.Q:pct_esp_vert_2011ct
## 0.0017542437
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_mat_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_gr
Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.062}\) = 0.94 (95%CI: 0.913 – 0.968)
Coefficient interpretation (marginal + interaction) as relative risk.
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(wSCOREMAT.2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_2011ct'] + i_e["wSCOREMAT.2011.Q"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_2011ct = quantile(clean_bei$df$pct_esp_vert_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['wSCOREMAT.2011.Q'] + i_e["pct_esp_vert_2011ct"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.085884 0.0497948 -41.890 0.000e+00
## wSCOREMAT.2011.Q -0.119070 0.0160742 -7.408 1.287e-13
## pct_esp_vert_high_2011ct 0.028240 0.0016388 17.232 0.000e+00
## wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct 0.005592 0.0006923 8.077 6.661e-16
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1379013
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.0105
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1791.84
ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.2011.Q")
## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q
## -0.15066242 -0.08758274
ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_high_2011ct")
## lower pct_esp_vert_high_2011ct upper pct_esp_vert_high_2011ct
## 0.02502735 0.03147717
ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct")
## lower wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct
## 0.004235836
## upper wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct
## 0.006962135
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_mat_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_tc
Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.119}\) = 0.888 (95%CI: 0.86 – 0.916)
Coefficient interpretation (marginal + interaction) as relative risk.
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(wSCOREMAT.2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_high_2011ct'] + i_e["wSCOREMAT.2011.Q"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_high_2011ct = quantile(clean_bei$df$pct_esp_vert_high_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['wSCOREMAT.2011.Q'] + i_e["pct_esp_vert_high_2011ct"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -3.2206 0.9043 -3.5614
## vis_minority_2011.Q -0.4739 0.2291 -2.0684
## Bike_lane.by.street.2011ct 4.9368 4.1706 1.1837
## vis_minority_2011.Q:Bike_lane.by.street.2011ct -0.5239 0.9395 -0.5576
## p-value
## (Intercept) 0.0003689
## vis_minority_2011.Q 0.0386031
## Bike_lane.by.street.2011ct 0.2365265
## vis_minority_2011.Q:Bike_lane.by.street.2011ct 0.5770925
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1379043
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 6.547
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -145.4209
# p-value not far from 0.05, use LRT to refine test
res.lmm.0 <- update(res.lmm.zero, . ~ . - vis_minority_2011.Q)
LRT(res.lmm.0, res.lmm.zero)
## chi2_LR df p_value
## p_v 3.522643 1 0.06053585
# NB: not converging
# ci.lmm.0.1 <- confint(res.lmm.0, "vis_minority_2011.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Bike_lane.by.street.2011ct")
# ci.lmm.0.3 <- confint(res.lmm.0, "vis_minority_2011.Q:Bike_lane.by.street.2011ct")
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
plot(simout)
}
# Plotting interaction || MEANINGLESS
# .t <- clean_bei$df
# .t["predict"] <- as.vector(predict(res.lmm.0))
#
# ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
# geom_smooth(method = "lm", alpha=.1) +
# labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -2.635e+00 0.086877 -30.3275
## vis_minority_2011.Q -2.080e-02 0.024749 -0.8404
## Bike_lane.by.street.2011ct 5.415e-02 0.004618 11.7247
## vis_minority_2011.Q:Bike_lane.by.street.2011ct -2.446e-05 0.001491 -0.0164
## p-value
## (Intercept) 0.0000
## vis_minority_2011.Q 0.4007
## Bike_lane.by.street.2011ct 0.0000
## vis_minority_2011.Q:Bike_lane.by.street.2011ct 0.9869
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.150489
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1121
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1880.959
# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "vis_minority_2011.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Bike_lane.by.street.2011ct")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "vis_minority_2011.Q:Bike_lane.by.street.2011ct")
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))
g_vis_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
geom_smooth(method = "lm", alpha=.1, linetype = "dashed") +
labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_bk.hrdl
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4064.759
Interpretation of the vis_minority_2011.Q coefficient | Zero model (not significant): we have \(e^{-0.474}\) = 0.623 Interpretation of the wSCOREMAT.2011.Q coefficient | Count model (not significant): we have \(e^{-0.021}\) = 0.979
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.464947 0.0488681 -29.978 0.000e+00
## vis_minority_2011.Q -0.058931 0.0147090 -4.006 6.163e-05
## pct_esp_vert_2011ct 0.016177 0.0008521 18.984 0.000e+00
## vis_minority_2011.Q:pct_esp_vert_2011ct 0.001079 0.0002848 3.788 1.518e-04
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.136093
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004657
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2009.835
ci.lmm.1 <- confint(res.lmm, "vis_minority_2011.Q")
## lower vis_minority_2011.Q upper vis_minority_2011.Q
## -0.08777005 -0.02998448
ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_2011ct")
## lower pct_esp_vert_2011ct upper pct_esp_vert_2011ct
## 0.01449601 0.01788267
ci.lmm.3 <- confint(res.lmm, "vis_minority_2011.Q:pct_esp_vert_2011ct")
## lower vis_minority_2011.Q:pct_esp_vert_2011ct
## 0.0005204252
## upper vis_minority_2011.Q:pct_esp_vert_2011ct
## 0.0016392402
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_vis_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_gr
Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.059}\) = 0.943 (95%CI: 0.916 – 0.97)
Coefficient interpretation (marginal + interaction) as relative risk.
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(vis_minority_2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_2011ct'] + i_e["vis_minority_2011.Q"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_2011ct = quantile(clean_bei$df$pct_esp_vert_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['vis_minority_2011.Q'] + i_e["pct_esp_vert_2011ct"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -2.163285 0.0535686 -40.383
## vis_minority_2011.Q -0.076677 0.0166197 -4.614
## pct_esp_vert_high_2011ct 0.031788 0.0017899 17.759
## vis_minority_2011.Q:pct_esp_vert_high_2011ct 0.003336 0.0006397 5.215
## p-value
## (Intercept) 0.000e+00
## vis_minority_2011.Q 3.957e-06
## pct_esp_vert_high_2011ct 0.000e+00
## vis_minority_2011.Q:pct_esp_vert_high_2011ct 1.842e-07
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1379765
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01184
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1808.484
ci.lmm.1 <- confint(res.lmm, "vis_minority_2011.Q")
## lower vis_minority_2011.Q upper vis_minority_2011.Q
## -0.10934259 -0.04405318
ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_high_2011ct")
## lower pct_esp_vert_high_2011ct upper pct_esp_vert_high_2011ct
## 0.02828030 0.03533605
ci.lmm.3 <- confint(res.lmm, "vis_minority_2011.Q:pct_esp_vert_high_2011ct")
## lower vis_minority_2011.Q:pct_esp_vert_high_2011ct
## 0.002081562
## upper vis_minority_2011.Q:pct_esp_vert_high_2011ct
## 0.004598848
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_vis_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_tc
Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.077}\) = 0.926 (95%CI: 0.896 – 0.957)
Coefficient interpretation (marginal + interaction) as relative risk.
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(vis_minority_2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_high_2011ct'] + i_e["vis_minority_2011.Q"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_high_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_high_2011ct = quantile(clean_bei$df$pct_esp_vert_high_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['vis_minority_2011.Q'] + i_e["pct_esp_vert_high_2011ct"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_high_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e
Gentrified CT between 2011 and 2016
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -4.9123 0.5013 -9.7998
## gentrified_2016_2011TRUE 0.1917 0.5944 0.3226
## Bike_lane.by.street.2011ct 3.4688 1.4200 2.4427
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct -1.6031 1.6952 -0.9457
## p-value
## (Intercept) 0.00000
## gentrified_2016_2011TRUE 0.74703
## Bike_lane.by.street.2011ct 0.01458
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct 0.34432
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1379976
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 6.356
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -146.8066
# NB: not converging
# ci.lmm.0.1 <- confint(res.lmm.0, "gentrified_2016_2011TRUE")
# ci.lmm.0.2 <- confint(res.lmm.0, "Bike_lane.by.street.2011ct")
# ci.lmm.0.3 <- confint(res.lmm.0, "gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct")
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
plot(simout)
}
# Plotting interaction || MEANINGLESS
# .t <- clean_bei$df
# .t["predict"] <- as.vector(predict(res.lmm.0))
#
# ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) +
# geom_smooth(method = "lm", alpha=.1) +
# labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -2.75518 0.046857 -58.800
## gentrified_2016_2011TRUE 0.20178 0.071981 2.803
## Bike_lane.by.street.2011ct 0.05753 0.002441 23.564
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct -0.01157 0.004373 -2.646
## p-value
## (Intercept) 0.000000
## gentrified_2016_2011TRUE 0.005060
## Bike_lane.by.street.2011ct 0.000000
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct 0.008136
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1504352
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1088
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1877.541
# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "gentrified_2016_2011TRUE")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Bike_lane.by.street.2011ct")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct")
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))
g_gen_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_bk.hrdl
# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4060.696
Interpretation of the gentrified_2016_2011 coefficient | Zero model (not significant): we have \(e^{0.192}\) = 1.211 Interpretation of the gentrified_2016_2011 coefficient | Count model (not significant): we have \(e^{0.202}\) = 1.224
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ gentrified_2016_2011 * pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ gentrified_2016_2011 * pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ gentrified_2016_2011 * pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -1.613014 0.0257795 -62.570
## gentrified_2016_2011TRUE -0.202554 0.0643593 -3.147
## pct_esp_vert_2011ct 0.018742 0.0004559 41.105
## gentrified_2016_2011TRUE:pct_esp_vert_2011ct 0.005263 0.0016871 3.120
## p-value
## (Intercept) 0.000000
## gentrified_2016_2011TRUE 0.001648
## pct_esp_vert_2011ct 0.000000
## gentrified_2016_2011TRUE:pct_esp_vert_2011ct 0.001811
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1349448
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004838
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2012.743
ci.lmm.1 <- confint(res.lmm, "gentrified_2016_2011TRUE")
## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE
## -0.32924546 -0.07617365
ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_2011ct")
## lower pct_esp_vert_2011ct upper pct_esp_vert_2011ct
## 0.01782698 0.01968221
ci.lmm.3 <- confint(res.lmm, "gentrified_2016_2011TRUE:pct_esp_vert_2011ct")
## lower gentrified_2016_2011TRUE:pct_esp_vert_2011ct
## 0.001943885
## upper gentrified_2016_2011TRUE:pct_esp_vert_2011ct
## 0.008578808
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_gen_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_gr
Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.203}\) = 0.817 (95%CI: 0.719 – 0.927)
Coefficient interpretation (marginal + interaction) as relative risk.
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(gentrified_2016_2011=0:1)
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_2011ct'] + i_e["gentrified_2016_2011"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_2011ct = quantile(clean_bei$df$pct_esp_vert_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['gentrified_2016_2011TRUE'] + i_e["pct_esp_vert_2011ct"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ gentrified_2016_2011 * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ gentrified_2016_2011 * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ gentrified_2016_2011 * pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -2.35730 0.0279340 -84.388
## gentrified_2016_2011TRUE -0.20459 0.0715883 -2.858
## pct_esp_vert_high_2011ct 0.03939 0.0009182 42.893
## gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct 0.01322 0.0035508 3.722
## p-value
## (Intercept) 0.0000000
## gentrified_2016_2011TRUE 0.0042654
## pct_esp_vert_high_2011ct 0.0000000
## gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct 0.0001974
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380802
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01222
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1814.314
ci.lmm.1 <- confint(res.lmm, "gentrified_2016_2011TRUE")
## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE
## -0.34539306 -0.06456076
ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_high_2011ct")
## lower pct_esp_vert_high_2011ct upper pct_esp_vert_high_2011ct
## 0.03756862 0.04126454
ci.lmm.3 <- confint(res.lmm, "gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct")
## lower gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct
## 0.006232907
## upper gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct
## 0.020169474
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_gen_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_tc
Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.205}\) = 0.815 (95%CI: 0.708 – 0.937)
Coefficient interpretation (marginal + interaction) as relative risk.
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(gentrified_2016_2011=0:1)
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_high_2011ct'] + i_e["gentrified_2016_2011"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e
# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_high_2011ct = quantile(clean_bei$df$pct_esp_vert_high_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['gentrified_2016_2011TRUE'] + i_e["pct_esp_vert_high_2011ct"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e
Following Yan and Behzad’s request, the exact same models as above are run, this time using SES quintiles as categories instead of ordinal values.
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ factor(wSCOREMAT.2011.Q) *
## Bike_lane.by.street.2011ct + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -6.438 9.722e-01 -6.6221
## factor(wSCOREMAT.2011.Q)2 1.926 1.005e+00 1.9166
## factor(wSCOREMAT.2011.Q)3 1.830 1.084e+00 1.6888
## factor(wSCOREMAT.2011.Q)4 1.713 1.083e+00 1.5818
## factor(wSCOREMAT.2011.Q)5 1.119 1.025e+00 1.0914
## Bike_lane.by.street.2011ct 54.449 1.422e+02 0.3830
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct -52.827 1.422e+02 -0.3716
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct 6167.495 3.011e+04 0.2048
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct 170.855 9.721e+02 0.1758
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct -52.766 1.422e+02 -0.3712
## p-value
## (Intercept) 3.541e-11
## factor(wSCOREMAT.2011.Q)2 5.529e-02
## factor(wSCOREMAT.2011.Q)3 9.126e-02
## factor(wSCOREMAT.2011.Q)4 1.137e-01
## factor(wSCOREMAT.2011.Q)5 2.751e-01
## Bike_lane.by.street.2011ct 7.017e-01
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct 7.102e-01
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct 8.377e-01
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct 8.605e-01
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct 7.105e-01
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.137743
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 5.117
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -142.7656
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
plot(simout)
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE
## (Intercept) -2.596117 0.077232
## factor(wSCOREMAT.2011.Q)2 -0.082531 0.095227
## factor(wSCOREMAT.2011.Q)3 -0.344225 0.101073
## factor(wSCOREMAT.2011.Q)4 -0.169518 0.101080
## factor(wSCOREMAT.2011.Q)5 0.019403 0.103981
## Bike_lane.by.street.2011ct 0.049586 0.003791
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct 0.002008 0.005510
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct 0.023129 0.006319
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct 0.006693 0.005783
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct -0.004647 0.006463
## t-value p-value
## (Intercept) -33.6146 0.0000000
## factor(wSCOREMAT.2011.Q)2 -0.8667 0.3861205
## factor(wSCOREMAT.2011.Q)3 -3.4057 0.0006599
## factor(wSCOREMAT.2011.Q)4 -1.6771 0.0935277
## factor(wSCOREMAT.2011.Q)5 0.1866 0.8519723
## Bike_lane.by.street.2011ct 13.0796 0.0000000
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct 0.3645 0.7154798
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct 3.6605 0.0002517
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct 1.1573 0.2471496
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct -0.7190 0.4721634
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1503961
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1058
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1874.29
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))
g_mat_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
geom_smooth(method = "lm", alpha=.1, linetype="dashed") +
labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_bk.hrdl
# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4058.111
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q)* pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q) * pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -1.4971348 0.0460467 -32.5134
## factor(wSCOREMAT.2011.Q)2 -0.1162630 0.0618502 -1.8797
## factor(wSCOREMAT.2011.Q)3 -0.0825597 0.0634354 -1.3015
## factor(wSCOREMAT.2011.Q)4 -0.1729263 0.0685627 -2.5222
## factor(wSCOREMAT.2011.Q)5 -0.4873193 0.0840536 -5.7977
## pct_esp_vert_2011ct 0.0171251 0.0007596 22.5460
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_2011ct 0.0014829 0.0010502 1.4119
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_2011ct 0.0009409 0.0011150 0.8438
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_2011ct 0.0029923 0.0013692 2.1854
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_2011ct 0.0119756 0.0022691 5.2778
## p-value
## (Intercept) 0.000e+00
## factor(wSCOREMAT.2011.Q)2 6.014e-02
## factor(wSCOREMAT.2011.Q)3 1.931e-01
## factor(wSCOREMAT.2011.Q)4 1.166e-02
## factor(wSCOREMAT.2011.Q)5 6.722e-09
## pct_esp_vert_2011ct 0.000e+00
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_2011ct 1.580e-01
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_2011ct 3.988e-01
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_2011ct 2.886e-02
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_2011ct 1.308e-07
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1367616
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004266
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2000.254
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_mat_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_gr
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) * pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -2.195912 0.043413 -50.582
## factor(wSCOREMAT.2011.Q)2 -0.170865 0.062582 -2.730
## factor(wSCOREMAT.2011.Q)3 -0.176232 0.065792 -2.679
## factor(wSCOREMAT.2011.Q)4 -0.331512 0.070623 -4.694
## factor(wSCOREMAT.2011.Q)5 -0.749906 0.089324 -8.395
## pct_esp_vert_high_2011ct 0.033787 0.001234 27.370
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_high_2011ct 0.007058 0.002075 3.401
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_high_2011ct 0.006521 0.002386 2.734
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_high_2011ct 0.015541 0.002997 5.185
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_high_2011ct 0.042632 0.005363 7.949
## p-value
## (Intercept) 0.000e+00
## factor(wSCOREMAT.2011.Q)2 6.328e-03
## factor(wSCOREMAT.2011.Q)3 7.393e-03
## factor(wSCOREMAT.2011.Q)4 2.677e-06
## factor(wSCOREMAT.2011.Q)5 0.000e+00
## pct_esp_vert_high_2011ct 0.000e+00
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_high_2011ct 6.706e-04
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_high_2011ct 6.262e-03
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_high_2011ct 2.158e-07
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_high_2011ct 1.887e-15
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1379822
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.009176
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1780.152
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_mat_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_tc
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ factor(vis_minority_2011.Q) *
## Bike_lane.by.street.2011ct + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE
## (Intercept) -5.705 1.227
## factor(vis_minority_2011.Q)2 1.711 1.401
## factor(vis_minority_2011.Q)3 2.099 1.340
## factor(vis_minority_2011.Q)4 -0.692 1.451
## factor(vis_minority_2011.Q)5 -1.040 1.346
## Bike_lane.by.street.2011ct 63.221 369.060
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct 343.774 8009.562
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct -60.791 369.066
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct 418.608 4648.758
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct -60.552 369.061
## t-value p-value
## (Intercept) -4.65042 3.313e-06
## factor(vis_minority_2011.Q)2 1.22111 2.220e-01
## factor(vis_minority_2011.Q)3 1.56680 1.172e-01
## factor(vis_minority_2011.Q)4 -0.47704 6.333e-01
## factor(vis_minority_2011.Q)5 -0.77282 4.396e-01
## Bike_lane.by.street.2011ct 0.17130 8.640e-01
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct 0.04292 9.658e-01
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct -0.16472 8.692e-01
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct 0.09005 9.282e-01
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct -0.16407 8.697e-01
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1378827
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 11.6
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -137.7682
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
plot(simout)
}
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE
## (Intercept) -2.600489 0.078225
## factor(vis_minority_2011.Q)2 -0.114586 0.096774
## factor(vis_minority_2011.Q)3 -0.129058 0.100608
## factor(vis_minority_2011.Q)4 -0.250708 0.104275
## factor(vis_minority_2011.Q)5 -0.028887 0.107428
## Bike_lane.by.street.2011ct 0.051362 0.004280
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct 0.003173 0.005676
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct 0.006144 0.006687
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct 0.012759 0.006636
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct -0.005974 0.006801
## t-value p-value
## (Intercept) -33.2439 0.00000
## factor(vis_minority_2011.Q)2 -1.1840 0.23639
## factor(vis_minority_2011.Q)3 -1.2828 0.19957
## factor(vis_minority_2011.Q)4 -2.4043 0.01620
## factor(vis_minority_2011.Q)5 -0.2689 0.78801
## Bike_lane.by.street.2011ct 11.9997 0.00000
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct 0.5590 0.57613
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct 0.9189 0.35817
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct 1.9229 0.05449
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct -0.8785 0.37969
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1504866
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1102
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1876.529
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))
g_vis_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
geom_smooth(method = "lm", alpha=.1, linetype = "dashed") +
labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_bk.hrdl
# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4052.594
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) -1.501387 0.0461251 -32.550
## factor(vis_minority_2011.Q)2 -0.126758 0.0617600 -2.052
## factor(vis_minority_2011.Q)3 -0.119219 0.0673443 -1.770
## factor(vis_minority_2011.Q)4 -0.150541 0.0621047 -2.424
## factor(vis_minority_2011.Q)5 -0.327300 0.0697281 -4.694
## pct_esp_vert_2011ct 0.017003 0.0007648 22.233
## factor(vis_minority_2011.Q)2:pct_esp_vert_2011ct 0.001859 0.0010595 1.755
## factor(vis_minority_2011.Q)3:pct_esp_vert_2011ct 0.002070 0.0012848 1.611
## factor(vis_minority_2011.Q)4:pct_esp_vert_2011ct 0.002399 0.0011290 2.125
## factor(vis_minority_2011.Q)5:pct_esp_vert_2011ct 0.006502 0.0014860 4.375
## p-value
## (Intercept) 0.000e+00
## factor(vis_minority_2011.Q)2 4.013e-02
## factor(vis_minority_2011.Q)3 7.668e-02
## factor(vis_minority_2011.Q)4 1.535e-02
## factor(vis_minority_2011.Q)5 2.680e-06
## pct_esp_vert_2011ct 0.000e+00
## factor(vis_minority_2011.Q)2:pct_esp_vert_2011ct 7.932e-02
## factor(vis_minority_2011.Q)3:pct_esp_vert_2011ct 1.071e-01
## factor(vis_minority_2011.Q)4:pct_esp_vert_2011ct 3.357e-02
## factor(vis_minority_2011.Q)5:pct_esp_vert_2011ct 1.213e-05
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1358167
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.0045
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2006.416
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_vis_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_gr
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE
## (Intercept) -2.2494769 0.047193
## factor(vis_minority_2011.Q)2 -0.1218389 0.064577
## factor(vis_minority_2011.Q)3 0.0146893 0.066527
## factor(vis_minority_2011.Q)4 -0.2778787 0.065932
## factor(vis_minority_2011.Q)5 -0.3623408 0.077380
## pct_esp_vert_high_2011ct 0.0357059 0.001466
## factor(vis_minority_2011.Q)2:pct_esp_vert_high_2011ct 0.0041755 0.002176
## factor(vis_minority_2011.Q)3:pct_esp_vert_high_2011ct -0.0007537 0.002489
## factor(vis_minority_2011.Q)4:pct_esp_vert_high_2011ct 0.0108121 0.002458
## factor(vis_minority_2011.Q)5:pct_esp_vert_high_2011ct 0.0171154 0.003341
## t-value p-value
## (Intercept) -47.6658 0.000e+00
## factor(vis_minority_2011.Q)2 -1.8867 5.920e-02
## factor(vis_minority_2011.Q)3 0.2208 8.252e-01
## factor(vis_minority_2011.Q)4 -4.2147 2.502e-05
## factor(vis_minority_2011.Q)5 -4.6826 2.832e-06
## pct_esp_vert_high_2011ct 24.3551 0.000e+00
## factor(vis_minority_2011.Q)2:pct_esp_vert_high_2011ct 1.9187 5.502e-02
## factor(vis_minority_2011.Q)3:pct_esp_vert_high_2011ct -0.3028 7.620e-01
## factor(vis_minority_2011.Q)4:pct_esp_vert_high_2011ct 4.3986 1.089e-05
## factor(vis_minority_2011.Q)5:pct_esp_vert_high_2011ct 5.1230 3.007e-07
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380247
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01063
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1799.683
if (plot.resid) {
simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
plot(simout)
}
# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))
g_vis_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) +
geom_smooth(method = "lm", alpha=.1) +
labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_tc
No new models here, as gentrified variable is already used as a pure categorical variable.
As discussed with Yan and Behzad, we rerun the UI models, without the interaction term \(SES \times UC_{2011}\)
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q + Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ wSCOREMAT.2011.Q + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ wSCOREMAT.2011.Q + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -4.939506 0.8771 -5.63191 1.782e-08
## wSCOREMAT.2011.Q 0.003004 0.2117 0.01419 9.887e-01
## Bike_lane.by.street.2011ct 2.499779 0.6984 3.57955 3.442e-04
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1377147
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 5.786
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -149.7332
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.677404 0.068370 -39.1606 0.0000
## wSCOREMAT.2011.Q -0.007146 0.016550 -0.4318 0.6659
## Bike_lane.by.street.2011ct 0.054017 0.002135 25.2953 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1502893
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1125
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1884.431
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4078.329
Transforming the model coefficients (continuous):
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.007 |
| wSCOREMAT.2011.Q | 1.003 |
| Bike_lane.by.street.2011ct | 12.180 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.069 |
| wSCOREMAT.2011.Q | 0.993 |
| Bike_lane.by.street.2011ct | 1.056 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ wSCOREMAT.2011.Q + pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ wSCOREMAT.2011.Q + pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ wSCOREMAT.2011.Q + pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.58398 0.0342580 -46.24 0.0000
## wSCOREMAT.2011.Q -0.01166 0.0058884 -1.98 0.0477
## pct_esp_vert_2011ct 0.01875 0.0004615 40.63 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1359291
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004773
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2016.65
# Wald-test p-value for wSCOREMAT.2011.Q close to threshold, using LRT to improve accuracy
res.lmm.0 <- update(res.lmm, . ~ . -wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm)
## chi2_LR df p_value
## p_v 3.792348 1 0.05148739
Transforming the model coefficients (continuous):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.205 |
| wSCOREMAT.2011.Q | 0.988 |
| pct_esp_vert_2011ct | 1.019 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q + pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q + pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q + pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.32838 0.0424926 -54.795 0.0000
## wSCOREMAT.2011.Q -0.01022 0.0090720 -1.126 0.2601
## pct_esp_vert_high_2011ct 0.03961 0.0009664 40.985 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1378526
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01294
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1823.392
Transforming the model coefficients (continuous):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.097 |
| wSCOREMAT.2011.Q | 0.990 |
| pct_esp_vert_high_2011ct | 1.040 |
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q + Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ vis_minority_2011.Q + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ vis_minority_2011.Q + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -3.072 0.9195 -3.341 0.0008351
## vis_minority_2011.Q -0.516 0.2337 -2.208 0.0272372
## Bike_lane.by.street.2011ct 2.861 0.8191 3.493 0.0004782
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1378998
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 6.963
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -145.6995
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.63391 0.069833 -37.717 0.0000
## vis_minority_2011.Q -0.02109 0.017257 -1.222 0.2216
## Bike_lane.by.street.2011ct 0.05408 0.002161 25.023 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1504889
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1121
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1880.959
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4063.317
Transforming the model coefficients (continuous):
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.046 |
| vis_minority_2011.Q | 0.597 |
| Bike_lane.by.street.2011ct | 17.474 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.072 |
| vis_minority_2011.Q | 0.979 |
| Bike_lane.by.street.2011ct | 1.056 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ vis_minority_2011.Q + pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ vis_minority_2011.Q + pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ vis_minority_2011.Q + pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.608537 0.0323275 -49.757 0.0000
## vis_minority_2011.Q -0.007267 0.0057005 -1.275 0.2023
## pct_esp_vert_2011ct 0.018988 0.0004481 42.376 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1358025
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004891
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2016.956
Transforming the model coefficients (continuous):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.200 |
| vis_minority_2011.Q | 0.993 |
| pct_esp_vert_2011ct | 1.019 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ vis_minority_2011.Q + pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ vis_minority_2011.Q + pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ vis_minority_2011.Q + pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.356969 0.0402505 -58.5575 0.0000
## vis_minority_2011.Q -0.002609 0.0088099 -0.2961 0.7672
## pct_esp_vert_high_2011ct 0.039985 0.0009402 42.5293 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.138067
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01303
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1821.833
Transforming the model coefficients (continuous):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.095 |
| vis_minority_2011.Q | 0.997 |
| pct_esp_vert_high_2011ct | 1.041 |
Gentrified CT between 2011 and 2016
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 + Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ gentrified_2016_2011 + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ gentrified_2016_2011 + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -4.82263 0.5165 -9.3369 0.0000000
## gentrified_2016_2011TRUE -0.07812 0.6074 -0.1286 0.8976627
## Bike_lane.by.street.2011ct 2.81386 0.8174 3.4424 0.0005766
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380204
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 7.002
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -147.7386
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.71508 0.04457 -60.912 0.0000
## gentrified_2016_2011TRUE 0.05633 0.04724 1.192 0.2331
## Bike_lane.by.street.2011ct 0.05436 0.00214 25.397 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.150439
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1116
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1881
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4067.478
Transforming the model coefficients:
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.008 |
| gentrified_2016_2011TRUE | 0.925 |
| Bike_lane.by.street.2011ct | 16.674 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.066 |
| gentrified_2016_2011TRUE | 1.058 |
| Bike_lane.by.street.2011ct | 1.056 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ gentrified_2016_2011 + pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ gentrified_2016_2011 + pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ gentrified_2016_2011 + pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.63121 0.0255404 -63.8681 0.000
## gentrified_2016_2011TRUE -0.01288 0.0206032 -0.6249 0.532
## pct_esp_vert_2011ct 0.01909 0.0004486 42.5552 0.000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1360243
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004979
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2017.555
Transforming the model coefficients:
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.196 |
| gentrified_2016_2011TRUE | 0.987 |
| pct_esp_vert_2011ct | 1.019 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ gentrified_2016_2011 + pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ gentrified_2016_2011 + pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ gentrified_2016_2011 + pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.37670 0.0279253 -85.109 0.0000
## gentrified_2016_2011TRUE 0.03580 0.0293925 1.218 0.2233
## pct_esp_vert_high_2011ct 0.04024 0.0009101 44.218 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380247
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01293
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1821.142
Transforming the model coefficients:
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.093 |
| gentrified_2016_2011TRUE | 1.036 |
| pct_esp_vert_high_2011ct | 1.041 |
As discussed with Yan and Behzad, we rerun the UI models, without the interaction term \(SES \times UC_{2011}\), SES as categories instaed of continuous
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) + Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ factor(wSCOREMAT.2011.Q) + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ factor(wSCOREMAT.2011.Q) +
## Bike_lane.by.street.2011ct + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -6.3848 1.0783 -5.9213 3.195e-09
## factor(wSCOREMAT.2011.Q)2 1.7251 1.1184 1.5425 1.230e-01
## factor(wSCOREMAT.2011.Q)3 2.1681 1.1958 1.8131 6.981e-02
## factor(wSCOREMAT.2011.Q)4 1.9986 1.1962 1.6709 9.475e-02
## factor(wSCOREMAT.2011.Q)5 0.7114 1.1433 0.6222 5.338e-01
## Bike_lane.by.street.2011ct 2.7278 0.8021 3.4007 6.722e-04
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1376078
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 7.864
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -147.1698
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.661013 0.062493 -42.58081 0.0000
## factor(wSCOREMAT.2011.Q)2 -0.039221 0.061482 -0.63793 0.5235
## factor(wSCOREMAT.2011.Q)3 -0.072592 0.063707 -1.13947 0.2545
## factor(wSCOREMAT.2011.Q)4 -0.074928 0.065152 -1.15005 0.2501
## factor(wSCOREMAT.2011.Q)5 -0.004295 0.073100 -0.05875 0.9531
## Bike_lane.by.street.2011ct 0.054073 0.002142 25.24000 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1503103
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.112
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1883.302
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4076.943
Transforming the model coefficients (categorical):
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.002 |
| factor(wSCOREMAT.2011.Q)2 | 5.613 |
| factor(wSCOREMAT.2011.Q)3 | 8.742 |
| factor(wSCOREMAT.2011.Q)4 | 7.379 |
| factor(wSCOREMAT.2011.Q)5 | 2.037 |
| Bike_lane.by.street.2011ct | 15.300 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.070 |
| factor(wSCOREMAT.2011.Q)2 | 0.962 |
| factor(wSCOREMAT.2011.Q)3 | 0.930 |
| factor(wSCOREMAT.2011.Q)4 | 0.928 |
| factor(wSCOREMAT.2011.Q)5 | 0.996 |
| Bike_lane.by.street.2011ct | 1.056 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q) + pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q) + pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q) + pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.58982 0.0305516 -52.037 0.000000
## factor(wSCOREMAT.2011.Q)2 -0.03055 0.0189787 -1.609 0.107508
## factor(wSCOREMAT.2011.Q)3 -0.02366 0.0203601 -1.162 0.245144
## factor(wSCOREMAT.2011.Q)4 -0.02242 0.0222606 -1.007 0.313861
## factor(wSCOREMAT.2011.Q)5 -0.06958 0.0269580 -2.581 0.009854
## pct_esp_vert_2011ct 0.01872 0.0004635 40.389 0.000000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1358061
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004722
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2014.941
Transforming the model coefficients (categorical):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.204 |
| factor(wSCOREMAT.2011.Q)2 | 0.970 |
| factor(wSCOREMAT.2011.Q)3 | 0.977 |
| factor(wSCOREMAT.2011.Q)4 | 0.978 |
| factor(wSCOREMAT.2011.Q)5 | 0.933 |
| pct_esp_vert_2011ct | 1.019 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) + pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) + pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) + pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.36271 0.0374207 -63.1390 0.00000
## factor(wSCOREMAT.2011.Q)2 0.03639 0.0293902 1.2383 0.21560
## factor(wSCOREMAT.2011.Q)3 0.02160 0.0318033 0.6791 0.49710
## factor(wSCOREMAT.2011.Q)4 0.02420 0.0345149 0.7011 0.48325
## factor(wSCOREMAT.2011.Q)5 -0.07695 0.0416488 -1.8475 0.06467
## pct_esp_vert_high_2011ct 0.03951 0.0009565 41.3107 0.00000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1378409
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.01245
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1819.267
Transforming the model coefficients (categorical):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.094 |
| factor(wSCOREMAT.2011.Q)2 | 1.037 |
| factor(wSCOREMAT.2011.Q)3 | 1.022 |
| factor(wSCOREMAT.2011.Q)4 | 1.024 |
| factor(wSCOREMAT.2011.Q)5 | 0.926 |
| pct_esp_vert_high_2011ct | 1.040 |
Bike lane length change
f <- Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) + Bike_lane.by.street.2011ct + offset(log(street_length.hm))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ factor(vis_minority_2011.Q) + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))
## formula: I(Bike_lane_total.hm.2016ct != 0) ~ factor(vis_minority_2011.Q) +
## Bike_lane.by.street.2011ct + offset(log(street_length.hm)) +
## adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -5.0763 1.144 -4.436 9.180e-06
## factor(vis_minority_2011.Q)2 1.4928 1.316 1.135 2.565e-01
## factor(vis_minority_2011.Q)3 1.3100 1.253 1.046 2.958e-01
## factor(vis_minority_2011.Q)4 0.1949 1.300 0.150 8.808e-01
## factor(vis_minority_2011.Q)5 -1.6542 1.276 -1.297 1.948e-01
## Bike_lane.by.street.2011ct 3.2327 1.019 3.174 1.505e-03
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1378309
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 11.06
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -143.4633
# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
filter(Bike_lane_total.hm.2016ct != 0) %>%
tidy_df(CT16, f)
res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) + Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))
## formula: Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) + Bike_lane.by.street.2011ct +
## offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.63738 0.061277 -43.0402 0.0000
## factor(vis_minority_2011.Q)2 -0.07088 0.059733 -1.1866 0.2354
## factor(vis_minority_2011.Q)3 -0.05550 0.063083 -0.8798 0.3790
## factor(vis_minority_2011.Q)4 -0.10233 0.066387 -1.5414 0.1232
## factor(vis_minority_2011.Q)5 -0.08137 0.074967 -1.0854 0.2778
## Bike_lane.by.street.2011ct 0.05428 0.002169 25.0282 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1504811
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.1117
## # of obs: 564; # of groups: ct_no, 564
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1880.374
# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)
# Combined AIC
"Combined Hurdle AIC:"
## [1] "Combined Hurdle AIC:"
-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)
## [1] 4063.675
Transforming the model coefficients (categorical):
# Zero model
round(exp(res.lmm.zero[["fixef"]]), 3) %>%
knitr::kable(caption="Zero model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.006 |
| factor(vis_minority_2011.Q)2 | 4.450 |
| factor(vis_minority_2011.Q)3 | 3.706 |
| factor(vis_minority_2011.Q)4 | 1.215 |
| factor(vis_minority_2011.Q)5 | 0.191 |
| Bike_lane.by.street.2011ct | 25.349 |
# Count model
round(exp(res.lmm.cnt[["fixef"]]), 3) %>%
knitr::kable(caption="Count model transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.072 |
| factor(vis_minority_2011.Q)2 | 0.932 |
| factor(vis_minority_2011.Q)3 | 0.946 |
| factor(vis_minority_2011.Q)4 | 0.903 |
| factor(vis_minority_2011.Q)5 | 0.922 |
| Bike_lane.by.street.2011ct | 1.056 |
Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) + pct_esp_vert_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) + pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) + pct_esp_vert_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -1.6151901 0.0296830 -54.41456 0.0000
## factor(vis_minority_2011.Q)2 -0.0190763 0.0195996 -0.97330 0.3304
## factor(vis_minority_2011.Q)3 -0.0005607 0.0216980 -0.02584 0.9794
## factor(vis_minority_2011.Q)4 -0.0183082 0.0215032 -0.85142 0.3945
## factor(vis_minority_2011.Q)5 -0.0396666 0.0254609 -1.55794 0.1192
## pct_esp_vert_2011ct 0.0189953 0.0004495 42.26227 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.135867
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.004871
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -2015.997
Transforming the model coefficients (categorical):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.199 |
| factor(vis_minority_2011.Q)2 | 0.981 |
| factor(vis_minority_2011.Q)3 | 0.999 |
| factor(vis_minority_2011.Q)4 | 0.982 |
| factor(vis_minority_2011.Q)5 | 0.961 |
| pct_esp_vert_2011ct | 1.019 |
Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level
f <- area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) + pct_esp_vert_high_2011ct + offset(log(ct_area_hct))
clean_bei <- tidy_df(bei_df_aoi, CT16, f)
res.lmm <- fitme(area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) + pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
family = poisson())
summary(res.lmm, details = c(p_value="Wald"))
## formula: area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) + pct_esp_vert_high_2011ct +
## offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value p-value
## (Intercept) -2.365558 0.0357953 -66.08573 0.0000
## factor(vis_minority_2011.Q)2 -0.001358 0.0302778 -0.04486 0.9642
## factor(vis_minority_2011.Q)3 0.033812 0.0330881 1.02189 0.3068
## factor(vis_minority_2011.Q)4 -0.018870 0.0337475 -0.55914 0.5761
## factor(vis_minority_2011.Q)5 -0.006233 0.0387782 -0.16073 0.8723
## pct_esp_vert_high_2011ct 0.039999 0.0009377 42.65560 0.0000
## --------------- Random effects ---------------
## Family: gaussian( link = identity )
## --- Correlation parameters:
## 1.rho
## 0.1380643
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## ct_no : 0.0129
## # of obs: 688; # of groups: ct_no, 688
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): -1820.503
Transforming the model coefficients (categorical):
round(exp(res.lmm[["fixef"]]), 3) %>%
knitr::kable(caption="Transformed coeffs")
| x | |
|---|---|
| (Intercept) | 0.094 |
| factor(vis_minority_2011.Q)2 | 0.999 |
| factor(vis_minority_2011.Q)3 | 1.034 |
| factor(vis_minority_2011.Q)4 | 0.981 |
| factor(vis_minority_2011.Q)5 | 0.994 |
| pct_esp_vert_high_2011ct | 1.041 |
No new models here, as gentrified variable is already used as a pure categorical variable.
NOTICE: This section needs to be updated
Presenting interaction terms (STROBE and more specifically here): best practice suggests we should include interaction terms for specific values of covariates.
| SES metric | BE metric | UC 2011 trend | UI 2011 -> 2016 trend (SES effect only, no interaction) |
|---|---|---|---|
| Pampalon / MAT 2011 (zero + cnt) | Bike lane | 0.790 (p=0.016) & 0.954 (p=0.031*) | 1.054 (p=0.807) & 0.992 (p=0.723) |
| Pampalon / MAT 2011 | Greenness | 0.911 (p<0.001) | 0.940 (p<0.001) |
| Pampalon / MAT 2011 | Tree canopy | 0.879 (p<0.001) | 0.888 (p<0.001) |
| Visible Minority 2011 (zero + cnt) | Bike lane | 0.809 (p=0.018) & 0.933 (p=0.016) | 0.623 (p=0.061*) & 0.979 (p=0.401) |
| Visible Minority 2011 | Greenness | 0.934 (p<0.001) | 0.943 (p<0.001) |
| Visible Minority 2011 | Tree canopy | 0.896 (p<0.001) | 0.926 (p<0.001) |
| Gentrified 2011-2016 (zero + cnt) | Bike lane | 0.936 (p=0.789) & 1.1 (p=0.236) | 1.211 (p=0.747) & 1.224 (p=0.005) |
| Gentrified 2011-2016 | Greenness | 0.822 (p<0.001) | 0.817 (p=0.002) |
| Gentrified 2011-2016 | Tree canopy | 0.86 (p=0.003) | 0.815 (p=0.004) |
p-values: Wald t-tests except when p-value close to threshold where Likelihood ratio test has been used (marked with *) Bold: significant
# Green space
g_bec_bike <- ggplot() +
geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane.by.street.2011ct)), lwd = 0) +
scale_fill_distiller(name = "Bike lane to\nstreet length (%)", palette = "YlGn", direction = 1) +
geom_sf(data=filter(CT16, interact_aoi), fill=NA, color="gray", size=.2, alpha=0.2) +
ggnewscale::new_scale_fill() +
geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=factor(units::drop_units(Bike_lane.by.street.2011ct) == 0)), lwd = 0.1) +
scale_fill_manual(name="No bike lane in CT", values = c("TRUE" = "white"), labels=NULL, na.translate = FALSE) +
labs(title = "Normalized bike lane length in 2011")
g_bei_bike <- ggplot() +
geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.by.street.2011.2016ct)), lwd=0) +
geom_sf(data=filter(CT16, interact_aoi), fill=NA, color="gray", size=.2, alpha=0.2) +
scale_fill_gradient2(name = "Change in normalized\nbike lane (%)")+
labs(title = "Bike lane change between 2011 and 2016")
# Green space
g_bec_green <- ggplot() +
geom_sf(data=inner_join(filter(CT16, interact_aoi), esp_vert_ct, by="GeoUID"), mapping = aes(fill=as.numeric(pct_esp_vert_2011)), lwd = 0) +
geom_sf(data=filter(CT16, interact_aoi), fill=NA, color="gray", size=.2, alpha=0.2) +
scale_fill_distiller(name = "Greenspace (%)", palette = "BrBG", direction = 1)+
labs(title = "Proportion of greenspace in 2011")
g_bei_green <- ggplot() +
geom_sf(data=inner_join(filter(CT16, interact_aoi), esp_vert_ct, by="GeoUID"), mapping = aes(fill=as.numeric(pct_esp_vert_diff_2011.2017)), lwd = 0) +
geom_sf(data=filter(CT16, interact_aoi), fill=NA, color="gray", size=.2, alpha=0.2) +
scale_fill_gradient2(name = "Greenspace\nchange (%)")+
labs(title = "Greenspace change between 2011 and 2017")
# Green space
g_bec_tree <- ggplot() +
geom_sf(data=inner_join(filter(CT16, interact_aoi), esp_vert_ct, by="GeoUID"), mapping = aes(fill=as.numeric(pct_esp_vert_high_2011)), lwd = 0) +
geom_sf(data=filter(CT16, interact_aoi), fill=NA, color="gray", size=.2, alpha=0.2) +
scale_fill_distiller(name = "Tree canopy (%)", palette = "BrBG", direction = 1)+
labs(title = "Proportion of tree canopy in 2011")
g_bei_tree <- ggplot() +
geom_sf(data=inner_join(filter(CT16, interact_aoi), esp_vert_ct, by="GeoUID"), mapping = aes(fill=as.numeric(pct_esp_vert_diff_high_2011.2017)), lwd = 0) +
geom_sf(data=filter(CT16, interact_aoi), fill=NA, color="gray", size=.2, alpha=0.2) +
scale_fill_gradient2(name = "Tree canopy\nchange (%)")+
labs(title = "Tree canopy change between 2011 and 2017")
# Create grid of maps
pcol1 <- plot_grid(
g_bec_bike + theme_void(),
g_bec_green + theme_void(),
g_bec_tree + theme_void(),
align = "vh",
labels = c("A", "B", "C"),
hjust = -1,
ncol = 1
)
pcol2 <- plot_grid(
g_bei_bike + theme_void(),
g_bei_green + theme_void(),
g_bei_tree + theme_void(),
align = "vh",
#labels = c("D", "E", "F"),
hjust = -1,
ncol = 1
)
# output all graphs
plot_grid(pcol1, pcol2, nrow=1)
# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
transmute(year = "2011",
zone = zone,
CT_UID = CT_UID,
Population = Population,
`Material score` = wSCOREMAT.2011,
`% visible minorities` = vis_minority_2011,
`gentrified` = gentrified_2016_2011,
`Bike lane` = Bike_lane_total.2011ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2011ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2011ct,
`Green space area (ha)` = area_esp_vert_2011ct,
`% green space within CT` = pct_esp_vert_2011ct,
`Tree canopy area (ha)` = area_esp_vert_high_2011ct,
`% tree canopy within CT` = pct_esp_vert_high_2011ct)
# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
transmute(year = "2016",
zone = zone,
CT_UID = CT_UID,
`Bike lane` = Bike_lane_total.2016ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2016ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2016ct,
`Green space area (ha)` = area_esp_vert_2017ct,
`% green space within CT` = pct_esp_vert_2017ct,
`Tree canopy area (ha)` = area_esp_vert_high_2017ct,
`% tree canopy within CT` = pct_esp_vert_high_2017ct)
.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year) %>%
summarize(across(where(is.numeric), median, na.rm=TRUE),
`% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
`Bike lane presence` = sum(`Bike lane`))
.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year) %>%
summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))
.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year) %>%
summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))
tab1 <- inner_join(.tab1_med, .tab1_p25, by="year", suffix = c("", " p25")) %>%
inner_join(.tab1_p75, by="year", suffix = c("", " p75")) %>%
transmute(year = year,
Population = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
`Material score` = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
`% visible minorities` = case_when(is.na(`% visible minorities`) ~ "",
TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
`% gentrified CTs` = case_when(year == "2016" ~ "",
TRUE ~ as.character(round(`% gentrified`, 1))),
`Bike lane presence (N)` = as.character(`Bike lane presence`),
`Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
`% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
`Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
`% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
`Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
`% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
tibble::column_to_rownames("year") %>%
t() %>%
as.data.frame()
tab1 %>% knitr::kable() %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| 2011 | 2016 | |
|---|---|---|
| Population | 3812 (2628 - 5157) | |
| Material score | 0.006 (-0.019 - 0.031) | |
| % visible minorities | 21.3 (13.9 - 34.1) | |
| % gentrified CTs | 29.5 | |
| Bike lane presence (N) | 503 | 592 |
| Bike lane length (hm) | 6.9 (0 - 18.5) | 10.7 (3 - 22.2) |
| % Bike lane to streets | 7 (0 - 14) | 10.9 (3.8 - 18.2) |
| Green space area (ha) | 23 (8 - 64) | 26 (10 - 68) |
| % green space within CT | 36 (24.7 - 47.4) | 40.7 (29.5 - 51) |
| Tree canopy area (ha) | 11 (4 - 24) | 12 (5 - 28) |
| % tree canopy within CT | 16.4 (11.8 - 22.6) | 20.2 (14.5 - 26.8) |
We look at Urban conditions in 2011 and 2016 within Q1 and Q5 of 2011 equity metrics
# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
filter(wSCOREMAT.2011.Q %in% c(1, 5) ) %>%
transmute(year = "2011",
quintile = case_when(wSCOREMAT.2011.Q == 1 ~ "Q1 Deprivation",
wSCOREMAT.2011.Q == 5 ~ "Q5 Deprivation"),
zone = zone,
CT_UID = CT_UID,
Population = Population,
`Material score` = wSCOREMAT.2011,
`% visible minorities` = vis_minority_2011,
`gentrified` = gentrified_2016_2011,
`Bike lane` = Bike_lane_total.2011ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2011ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2011ct,
`Green space area (ha)` = area_esp_vert_2011ct,
`% green space within CT` = pct_esp_vert_2011ct,
`Tree canopy area (ha)` = area_esp_vert_high_2011ct,
`% tree canopy within CT` = pct_esp_vert_high_2011ct)
# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
filter(wSCOREMAT.2011.Q %in% c(1, 5) ) %>%
transmute(year = "2016",
quintile = case_when(wSCOREMAT.2011.Q == 1 ~ "Q1 Deprivation",
wSCOREMAT.2011.Q == 5 ~ "Q5 Deprivation"),
zone = zone,
CT_UID = CT_UID,
`Bike lane` = Bike_lane_total.2016ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2016ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2016ct,
`Green space area (ha)` = area_esp_vert_2017ct,
`% green space within CT` = pct_esp_vert_2017ct,
`Tree canopy area (ha)` = area_esp_vert_high_2017ct,
`% tree canopy within CT` = pct_esp_vert_high_2017ct)
.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), median, na.rm=TRUE),
`% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
`Bike lane presence` = sum(`Bike lane`),
`Bike lane absence` = n() - sum(`Bike lane`))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
tab1 <- inner_join(.tab1_med, .tab1_p25, by=c("year", "quintile"), suffix = c("", " p25")) %>%
inner_join(.tab1_p75, by=c("year", "quintile"), suffix = c("", " p75")) %>%
transmute(year.quintile = paste(year, quintile),
Population = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
`Material score` = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
`% visible minorities` = case_when(is.na(`% visible minorities`) ~ "",
TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
`% gentrified CTs` = case_when(year == "2016" ~ "",
TRUE ~ as.character(round(`% gentrified`, 1))),
`Bike lane presence (N)` = as.character(`Bike lane presence`),
`Bike lane absence (N)` = as.character(`Bike lane absence`),
`Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
`% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
`Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
`% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
`Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
`% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
tibble::column_to_rownames("year.quintile") %>%
t() %>%
as.data.frame()
tab1 %>% knitr::kable() %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| 2011 Q1 Deprivation | 2011 Q5 Deprivation | 2016 Q1 Deprivation | 2016 Q5 Deprivation | |
|---|---|---|---|---|
| year | 2011 | 2011 | 2016 | 2016 |
| Population | 3149 (2026.8 - 5017.8) | 3904.5 (2978.8 - 5212.2) | ||
| Material score | -0.046 (-0.061 - -0.037) | 0.054 (0.045 - 0.07) | ||
| % visible minorities | 12.8 (7.9 - 21.1) | 37.1 (25.3 - 54.4) | ||
| % gentrified CTs | 13 | 40.6 | ||
| Bike lane presence (N) | 106 | 78 | 116 | 99 |
| Bike lane absence (N) | 32 | 60 | 22 | 39 |
| Bike lane length (hm) | 9.3 (0 - 31.9) | 2 (0 - 9.1) | 11.3 (3.3 - 40.3) | 5.6 (0 - 12.6) |
| % Bike lane to streets | 8.2 (0 - 18.2) | 2.4 (0 - 9.8) | 12.2 (4.5 - 21.7) | 8.1 (0 - 15.1) |
| Green space area (ha) | 31 (5.2 - 116) | 15.5 (7 - 24.8) | 36 (7 - 119.8) | 17 (8 - 28) |
| % green space within CT | 44.5 (26.1 - 54.6) | 29.1 (20.9 - 35.2) | 49.9 (31.2 - 59.1) | 33.3 (25 - 40.9) |
| Tree canopy area (ha) | 15 (4 - 48.8) | 7 (4 - 10) | 19.5 (5 - 54.5) | 8 (5 - 11.8) |
| % tree canopy within CT | 22 (14.5 - 34.2) | 12.8 (9 - 16) | 28.2 (18.2 - 37.4) | 15.9 (10.6 - 19.7) |
# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
filter(vis_minority_2011.Q %in% c(1, 5) ) %>%
transmute(year = "2011",
quintile = case_when(vis_minority_2011.Q == 1 ~ "Q1 Vis. Minorities",
vis_minority_2011.Q == 5 ~ "Q5 Vis. Minorities"),
zone = zone,
CT_UID = CT_UID,
Population = Population,
`Material score` = wSCOREMAT.2011,
`% visible minorities` = vis_minority_2011,
`gentrified` = gentrified_2016_2011,
`Bike lane` = Bike_lane_total.2011ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2011ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2011ct,
`Green space area (ha)` = area_esp_vert_2011ct,
`% green space within CT` = pct_esp_vert_2011ct,
`Tree canopy area (ha)` = area_esp_vert_high_2011ct,
`% tree canopy within CT` = pct_esp_vert_high_2011ct)
# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
filter(vis_minority_2011.Q %in% c(1, 5) ) %>%
transmute(year = "2016",
quintile = case_when(vis_minority_2011.Q == 1 ~ "Q1 Vis. Minorities",
vis_minority_2011.Q == 5 ~ "Q5 Vis. Minorities"),
zone = zone,
CT_UID = CT_UID,
`Bike lane` = Bike_lane_total.2016ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2016ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2016ct,
`Green space area (ha)` = area_esp_vert_2017ct,
`% green space within CT` = pct_esp_vert_2017ct,
`Tree canopy area (ha)` = area_esp_vert_high_2017ct,
`% tree canopy within CT` = pct_esp_vert_high_2017ct)
.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), median, na.rm=TRUE),
`% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
`Bike lane presence` = sum(`Bike lane`),
`Bike lane absence` = n() - sum(`Bike lane`))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
tab1 <- inner_join(.tab1_med, .tab1_p25, by=c("year", "quintile"), suffix = c("", " p25")) %>%
inner_join(.tab1_p75, by=c("year", "quintile"), suffix = c("", " p75")) %>%
transmute(year.quintile = paste(year, quintile),
Population = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
`Material score` = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
`% visible minorities` = case_when(is.na(`% visible minorities`) ~ "",
TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
`% gentrified CTs` = case_when(year == "2016" ~ "",
TRUE ~ as.character(round(`% gentrified`, 1))),
`Bike lane presence (N)` = as.character(`Bike lane presence`),
`Bike lane absence (N)` = as.character(`Bike lane absence`),
`Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
`% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
`Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
`% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
`Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
`% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
tibble::column_to_rownames("year.quintile") %>%
t() %>%
as.data.frame()
tab1 %>% knitr::kable() %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| 2011 Q1 Vis. Minorities | 2011 Q5 Vis. Minorities | 2016 Q1 Vis. Minorities | 2016 Q5 Vis. Minorities | |
|---|---|---|---|---|
| year | 2011 | 2011 | 2016 | 2016 |
| Population | 3101.5 (2295.8 - 4662.2) | 4139 (3391.2 - 5594.5) | ||
| Material score | -0.023 (-0.044 - 0.008) | 0.037 (0.011 - 0.064) | ||
| % visible minorities | 8.5 (6.2 - 10.6) | 48.6 (42.3 - 58) | ||
| % gentrified CTs | 29.7 | 29 | ||
| Bike lane presence (N) | 104 | 77 | 120 | 95 |
| Bike lane absence (N) | 34 | 61 | 18 | 43 |
| Bike lane length (hm) | 8.9 (0 - 30.9) | 2.2 (0 - 9.8) | 13.1 (5 - 35.9) | 5.2 (0 - 14.7) |
| % Bike lane to streets | 8.3 (0 - 15.4) | 2.2 (0 - 8.9) | 13.9 (7.4 - 21.2) | 6.8 (0 - 14.4) |
| Green space area (ha) | 26.5 (7 - 116.8) | 16 (7.2 - 32.5) | 29 (9 - 120.8) | 18.5 (9.2 - 36.8) |
| % green space within CT | 39.8 (26.5 - 56.7) | 30.4 (23 - 38.3) | 44.3 (33 - 60.1) | 35.8 (26.9 - 44.6) |
| Tree canopy area (ha) | 12.5 (4 - 48.5) | 7 (4 - 15) | 15 (5.2 - 53.8) | 9.5 (4 - 17) |
| % tree canopy within CT | 19.1 (14.6 - 29.8) | 14.5 (8.6 - 20.6) | 23.4 (18.1 - 34.1) | 17.4 (9.9 - 24.2) |
# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
transmute(year = "2011",
quintile = case_when(gentrified_2016_2011 ~ "Gentrified",
TRUE ~ "Not gentrified"),
zone = zone,
CT_UID = CT_UID,
Population = Population,
`Material score` = wSCOREMAT.2011,
`% visible minorities` = vis_minority_2011,
`gentrified` = gentrified_2016_2011,
`Bike lane` = Bike_lane_total.2011ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2011ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2011ct,
`Green space area (ha)` = area_esp_vert_2011ct,
`% green space within CT` = pct_esp_vert_2011ct,
`Tree canopy area (ha)` = area_esp_vert_high_2011ct,
`% tree canopy within CT` = pct_esp_vert_high_2011ct)
# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
transmute(year = "2016",
quintile = case_when(gentrified_2016_2011 ~ "Gentrified",
TRUE ~ "Not gentrified"),
zone = zone,
CT_UID = CT_UID,
`Bike lane` = Bike_lane_total.2016ct > 0,
`Bike lane length (hm)` = Bike_lane_total.2016ct/100,
`% Bike lane to streets` = Bike_lane.by.street.2016ct,
`Green space area (ha)` = area_esp_vert_2017ct,
`% green space within CT` = pct_esp_vert_2017ct,
`Tree canopy area (ha)` = area_esp_vert_high_2017ct,
`% tree canopy within CT` = pct_esp_vert_high_2017ct)
.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), median, na.rm=TRUE),
`% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
`Bike lane presence` = sum(`Bike lane`),
`Bike lane absence` = n() - sum(`Bike lane`))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
group_by(year, quintile) %>%
summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
tab1 <- inner_join(.tab1_med, .tab1_p25, by=c("year", "quintile"), suffix = c("", " p25")) %>%
inner_join(.tab1_p75, by=c("year", "quintile"), suffix = c("", " p75")) %>%
transmute(year.quintile = paste(year, quintile),
Population = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
`Material score` = case_when(year == "2016" ~ "",
TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
`% visible minorities` = case_when(is.na(`% visible minorities`) ~ "",
TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
`% gentrified CTs` = case_when(year == "2016" ~ "",
TRUE ~ as.character(round(`% gentrified`, 1))),
`Bike lane presence (N)` = as.character(`Bike lane presence`),
`Bike lane absence (N)` = as.character(`Bike lane absence`),
`Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
`% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
`Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
`% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
`Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
`% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
tibble::column_to_rownames("year.quintile") %>%
t() %>%
as.data.frame()
tab1 %>% knitr::kable() %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| 2011 Gentrified | 2011 Not gentrified | 2016 Gentrified | 2016 Not gentrified | |
|---|---|---|---|---|
| year | 2011 | 2011 | 2016 | 2016 |
| Population | 3159.5 (2266.8 - 4505.8) | 4043 (2862 - 5345) | ||
| Material score | 0.018 (-0.004 - 0.039) | 0.002 (-0.026 - 0.027) | ||
| % visible minorities | 21 (13.9 - 32.4) | 21.7 (14.5 - 34.7) | ||
| % gentrified CTs | 100 | 0 | ||
| Bike lane presence (N) | 130 | 373 | 171 | 421 |
| Bike lane absence (N) | 78 | 124 | 37 | 76 |
| Bike lane length (hm) | 3.7 (0 - 11) | 9.1 (0 - 23.5) | 8 (2.4 - 15.4) | 12.5 (3.7 - 26.8) |
| % Bike lane to streets | 5.9 (0 - 14.4) | 7.1 (0 - 13.7) | 13.3 (3.9 - 19.8) | 10.1 (3.7 - 17.6) |
| Green space area (ha) | 9 (4 - 22.2) | 34 (15 - 84) | 10.5 (5 - 25.2) | 38 (17 - 87) |
| % green space within CT | 27.5 (20.7 - 35.3) | 40.5 (29.1 - 51.3) | 32.8 (25.7 - 40.6) | 44.7 (34 - 54.6) |
| Tree canopy area (ha) | 5 (3 - 10) | 14 (6 - 33) | 6 (4 - 12) | 17 (8 - 37) |
| % tree canopy within CT | 15.2 (11.4 - 19.5) | 17 (11.9 - 24.2) | 19.4 (15 - 23.8) | 20.7 (14.1 - 28.4) |
# extract the legend from one of the plots
legend.row1 <- get_legend(
# create some space to the left of the legend
g_mat_gr + labs(colour='Material score\n(Quintiles)', fill='Material score\n(Quintiles)') + theme(legend.box.margin = margin(0, 0, 0, 12))
)
prow1.1 <- plot_grid(
g_mat_bk.hrdl + xlab("% Bike lane to streets 2011") + ylab("% Bike lane to streets 2016") + theme(legend.position = "none"),
g_mat_gr + xlab("% green space within CT 2011") + ylab("% green space within CT 2016") + theme(legend.position = "none"),
g_mat_tc + xlab("% tree canopy within CT 2011") + ylab("% tree canopy within CT 2016") + theme(legend.position = "none"),
align = "vh",
labels = c("A", "B", "C"),
hjust = -1,
nrow = 1
)
prow1 <- plot_grid(prow1.1, legend.row1, rel_widths = c(3, .4))
# extract the legend from one of the plots
legend.row2 <- get_legend(
# create some space to the left of the legend
g_vis_gr + labs(colour='Visible minorities\n(Quintiles)', fill='Visible minorities\n(Quintiles)') + theme(legend.box.margin = margin(0, 0, 0, 12))
)
prow2.1 <- plot_grid(
g_vis_bk.hrdl + xlab("% Bike lane to streets 2011") + ylab("% Bike lane to streets 2016") + theme(legend.position = "none"),
g_vis_gr + xlab("% green space within CT 2011") + ylab("% green space within CT 2016") + theme(legend.position = "none"),
g_vis_tc + xlab("% tree canopy within CT 2011") + ylab("% tree canopy within CT 2016") + theme(legend.position = "none"),
align = "vh",
labels = c("D", "E", "F"),
hjust = -1,
nrow = 1
)
prow2 <- plot_grid(prow2.1, legend.row2, rel_widths = c(3, .4))
# extract the legend from one of the plots
legend.row3 <- get_legend(
# create some space to the left of the legend
g_gen_gr + labs(colour='Gentrified', fill='Gentrified') + theme(legend.box.margin = margin(0, 0, 0, 12))
)
prow3.1 <- plot_grid(
g_gen_bk.hrdl + xlab("% Bike lane to streets 2011") + ylab("% Bike lane to streets 2016") + theme(legend.position = "none"),
g_gen_gr + xlab("% green space within CT 2011") + ylab("% green space within CT 2016") + theme(legend.position = "none"),
g_gen_tc + xlab("% tree canopy within CT 2011") + ylab("% tree canopy within CT 2016") + theme(legend.position = "none"),
align = "vh",
labels = c("G", "H", "I"),
hjust = -1,
nrow = 1
)
prow3 <- plot_grid(prow3.1, legend.row3, rel_widths = c(3, .4))
# output all graphs
plot_grid(prow1, prow2, prow3, ncol=1)
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cancensus_0.5.0 stringr_1.4.0 DHARMa_0.4.5 stargazer_5.2.2
## [5] spaMM_3.10.0 cachem_1.0.6 memoise_2.0.1 spatialreg_1.2-1
## [9] spdep_1.2-2 spData_2.0.1 sp_1.4-6 lme4_1.1-28
## [13] Matrix_1.4-0 cowplot_1.1.1 biscale_0.2.0 openxlsx_4.2.5
## [17] RPostgres_1.4.3 DBI_1.1.2 ggmap_3.0.0 ggplot2_3.3.5
## [21] stars_0.5-5 abind_1.4-5 sf_1.0-6 tidyr_1.2.0
## [25] dplyr_1.0.8
##
## loaded via a namespace (and not attached):
## [1] ggnewscale_0.4.7 minqa_1.2.4 colorspace_2.0-3
## [4] rjson_0.2.21 deldir_1.0-6 ellipsis_0.3.2
## [7] class_7.3-20 rstudioapi_0.13 proxy_0.4-26
## [10] farver_2.1.0 bit64_4.0.5 fansi_1.0.2
## [13] lubridate_1.8.0 xml2_1.3.3 codetools_0.2-18
## [16] splines_4.1.2 knitr_1.37 jsonlite_1.8.0
## [19] nloptr_2.0.0 png_0.1-7 compiler_4.1.2
## [22] httr_1.4.2 assertthat_0.2.1 fastmap_1.1.0
## [25] cli_3.2.0 s2_1.0.7 htmltools_0.5.2
## [28] tools_4.1.2 coda_0.19-4 gtable_0.3.0
## [31] glue_1.6.2 wk_0.6.0 gmodels_2.18.1
## [34] Rcpp_1.0.8 slam_0.1-50 jquerylib_0.1.4
## [37] raster_3.5-15 vctrs_0.3.8 svglite_2.1.0
## [40] gdata_2.18.0 nlme_3.1-155 lwgeom_0.2-8
## [43] xfun_0.29 rvest_1.0.2 lifecycle_1.0.1
## [46] gtools_3.9.2 terra_1.5-21 ROI_1.0-0
## [49] LearnBayes_2.15.1 MASS_7.3-55 scales_1.1.1
## [52] hms_1.1.1 parallel_4.1.2 expm_0.999-6
## [55] RColorBrewer_1.1-2 yaml_2.3.5 pbapply_1.5-0
## [58] sass_0.4.0 stringi_1.7.6 highr_0.9
## [61] e1071_1.7-9 boot_1.3-28 zip_2.2.0
## [64] systemfonts_1.0.4 RgoogleMaps_1.4.5.3 rlang_1.0.1
## [67] pkgconfig_2.0.3 bitops_1.0-7 evaluate_0.15
## [70] lattice_0.20-45 purrr_0.3.4 labeling_0.4.2
## [73] bit_4.0.4 tidyselect_1.1.2 geojsonsf_2.0.1
## [76] plyr_1.8.6 magrittr_2.0.2 R6_2.5.1
## [79] generics_0.1.2 mgcv_1.8-39 pillar_1.7.0
## [82] withr_2.4.3 units_0.8-0 tibble_3.1.6
## [85] crayon_1.5.0 KernSmooth_2.23-20 utf8_1.2.2
## [88] rmarkdown_2.11 jpeg_0.1-9 grid_4.1.2
## [91] blob_1.2.2 webshot_0.5.2 digest_0.6.29
## [94] classInt_0.4-3 numDeriv_2016.8-1.1 munsell_0.5.0
## [97] viridisLite_0.4.0 kableExtra_1.3.4 registry_0.5-1
## [100] bslib_0.3.1